LUCAS¶
LUCAS (Land Use/Cover Area Frame Statistical Survey) je datová sada, která poskytuje podrobné informace o využití půdy, půdním pokryvu a agroenvironmentálních proměnných v celé Evropské unii. Datová sada je doplňována prostřednictvím terénních šetření prováděných každé tři roky. Dataset je poskytován Eurostatem v režimu otevřených dat. Má zásadní význam pro monitorování životního prostředí a výzkum, zejména pro pochopení změn ve využívání půdy v průběhu času.
Projekt ST_LUCAS poskytuje standardizované rozhraní pro přístup k harmonizovaným a časoprostorově agregovaným datům LUCAS. Data jsou poskytována jako služba OGC Web Feature Service. Pro usnadnění přístupu k datům vznikl:
- QGIS Plugin (ST_LUCAS Download Manager)
- Python API
Nejprve si vykoušíme stažení dat pomocí QGIS pluginu. Nastavíme mapové okno na naše zájmové území a na základě něho stáhneme data LUCAS.
Projekt ST_LUCAS nabízí také Python API (viz návody). Balíček st_lucas si nejprve doinstalujeme (a posléze restartujeme kernel Jupyter notebooku).
In [ ]:
Copied!
!pip install st_lucas
!pip install st_lucas
Souřadnice minimálního ohraničujícího obdélníku získáme z prostředí QGISu (pozor na prohození souřadnic).
In [2]:
Copied!
from st_lucas import LucasRequest, LucasIO
request = LucasRequest()
request.bbox = (4587768, 3041298, 4610252, 3057649)
# request.years = [2018]
request.group = 'LC_LU'
lucasio = LucasIO()
lucasio.download(request)
lucasio.count()
from st_lucas import LucasRequest, LucasIO
request = LucasRequest()
request.bbox = (4587768, 3041298, 4610252, 3057649)
# request.years = [2018]
request.group = 'LC_LU'
lucasio = LucasIO()
lucasio.download(request)
lucasio.count()
2025-10-02 12:51:56,094 - LUCAS - INFO - io.download - Download process successfuly finished. Size of downloaded data: 318kb
Out[2]:
150
Objekt se staženými LUCAS body převedeme do GeoDataFrame pomocí metody to_geopandas():
In [3]:
Copied!
lucas_data = lucasio.to_geopandas()
lucas_data.head(5)
lucas_data = lucasio.to_geopandas()
lucas_data.head(5)
Out[3]:
| point_id | nuts0 | nuts1 | nuts2 | nuts3 | survey_date | car_latitude | car_longitude | car_ew | gps_proj | ... | photo_point | photo_north | photo_east | photo_south | photo_west | crop_residues | transect | ex_ante | survey_year | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 45943042 | CZ | CZ0 | CZ04 | None | 2006-05-16 | NaN | NaN | NaN | 1 | ... | NaN | NaN | NaN | NaN | NaN | NaN | None | NaN | 2006 | POINT (4594000.37 3042000.543) |
| 1 | 45963044 | CZ | CZ0 | CZ04 | None | 2006-05-15 | NaN | NaN | NaN | 1 | ... | NaN | NaN | NaN | NaN | NaN | NaN | None | NaN | 2006 | POINT (4596000.851 3043999.936) |
| 2 | 46003048 | CZ | CZ0 | CZ04 | None | 2006-05-15 | NaN | NaN | NaN | 1 | ... | NaN | NaN | NaN | NaN | NaN | NaN | None | NaN | 2006 | POINT (4600000.672 3048000.167) |
| 3 | 46103050 | CZ | CZ0 | CZ04 | None | 2006-05-15 | NaN | NaN | NaN | 1 | ... | NaN | NaN | NaN | NaN | NaN | NaN | None | NaN | 2006 | POINT (4610000.128 3050000.214) |
| 4 | 46103054 | CZ | CZ0 | CZ04 | None | 2006-05-15 | NaN | NaN | NaN | 1 | ... | NaN | NaN | NaN | NaN | NaN | NaN | None | NaN | 2006 | POINT (4610000.259 3053999.489) |
5 rows × 74 columns
In [4]:
Copied!
len(lucas_data.columns)
len(lucas_data.columns)
Out[4]:
74
In [5]:
Copied!
lucas_data.explore()
lucas_data.explore()
Out[5]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Spočítejme si jednoduchou statistiku opakovaného měření.
In [6]:
Copied!
lucas_data["survey_year"].unique()
lucas_data["survey_year"].unique()
Out[6]:
array([2006, 2009, 2012, 2015, 2018], dtype=int32)
In [7]:
Copied!
lucas_data.groupby("survey_year").count()["point_id"]
lucas_data.groupby("survey_year").count()["point_id"]
Out[7]:
survey_year 2006 33 2009 27 2012 29 2015 30 2018 31 Name: point_id, dtype: int64
Zaměříme se podrobněji na atribut lc1_h (Harmonized Land Cover 1 to 2018 nomenclature: nomenklatura)
In [8]:
Copied!
lucas_data.groupby("lc1_h").count()["point_id"].sort_values(ascending=False).head()
lucas_data.groupby("lc1_h").count()["point_id"].sort_values(ascending=False).head()
Out[8]:
lc1_h B11 33 B13 21 C10 17 E20 11 B32 11 Name: point_id, dtype: int64
In [9]:
Copied!
# pouze pro rok 2018
lucas_data[lucas_data["survey_year"] == 2018].groupby("lc1_h").\
count()["point_id"].sort_values(ascending=False).head()
# pouze pro rok 2018
lucas_data[lucas_data["survey_year"] == 2018].groupby("lc1_h").\
count()["point_id"].sort_values(ascending=False).head()
Out[9]:
lc1_h B11 9 B32 4 C10 4 B13 3 A21 2 Name: point_id, dtype: int64
Poznámka: ST_LUCAS umožňuje získat prostorově agregovaná data přímo:
In [10]:
Copied!
request.st_aggregated = True
request.years = [2012, 2015, 2018]
lucasio = LucasIO()
lucasio.download(request)
lucasio.to_geopandas()[["point_id", "lc1_h_2012", "lc1_h_2015", "lc1_h_2018"]].head(5)
request.st_aggregated = True
request.years = [2012, 2015, 2018]
lucasio = LucasIO()
lucasio.download(request)
lucasio.to_geopandas()[["point_id", "lc1_h_2012", "lc1_h_2015", "lc1_h_2018"]].head(5)
2025-10-02 12:51:58,334 - LUCAS - INFO - io.download - Download process successfuly finished. Size of downloaded data: 346kb
Out[10]:
| point_id | lc1_h_2012 | lc1_h_2015 | lc1_h_2018 | |
|---|---|---|---|---|
| 0 | 45963050 | B71 | B71 | None |
| 1 | 45963054 | C10 | C10 | None |
| 2 | 46063054 | C10 | C10 | None |
| 3 | 45943046 | B11 | B11 | B11 |
| 4 | 46003044 | C22 | C22 | C22 |
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.
V našem případě se zaměříme na produkt L2A (snímky, které byly korigovány na atmosférické vlivy).
Stažení dat¶
In [11]:
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/yusu/venv/lib/python3.13/site-packages/eodag/api/core.py:30: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81. import pkg_resources /home/martin/git/k155cvut/yusu/venv/lib/python3.13/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
2025-10-02 12:52:01,508 - eodag.config - INFO - config.override_config_from_file - Loading user configuration from: /home/martin/.config/eodag/eodag.yml 2025-10-02 12:52:01,518 - eodag.config - INFO - config.override_config_from_mapping - wekeo: unknown provider found in user conf, trying to use provided configuration 2025-10-02 12:52:01,519 - eodag.config - WARNING - config.override_config_from_mapping - wekeo skipped: could not be loaded from user configuration 2025-10-02 12:52:01,523 - eodag.core - INFO - core._prune_providers_list - aws_eos: provider needing auth for search has been pruned because no credentials could be found 2025-10-02 12:52:01,524 - eodag.core - INFO - core._prune_providers_list - meteoblue: provider needing auth for search has been pruned because no credentials could be found 2025-10-02 12:52:01,524 - eodag.core - INFO - core._prune_providers_list - hydroweb_next: provider needing auth for search has been pruned because no credentials could be found 2025-10-02 12:52:01,524 - eodag.core - INFO - core._prune_providers_list - wekeo_main: provider needing auth for search has been pruned because no credentials could be found 2025-10-02 12:52:01,525 - eodag.core - INFO - core._prune_providers_list - wekeo_ecmwf: provider needing auth for search has been pruned because no credentials could be found 2025-10-02 12:52:01,525 - eodag.core - INFO - core._prune_providers_list - wekeo_cmems: provider needing auth for search has been pruned because no credentials could be found 2025-10-02 12:52:01,525 - eodag.core - INFO - core._prune_providers_list - creodias_s3: provider needing auth for search has been pruned because no credentials could be found 2025-10-02 12:52:01,525 - eodag.core - INFO - core._prune_providers_list - dedl: provider needing auth for search has been pruned because no credentials could be found 2025-10-02 12:52:01,541 - eodag.core - INFO - core.set_locations_conf - Locations configuration loaded from /home/martin/.config/eodag/locations.yml
Out[11]:
['astraea_eod', 'cop_dataspace', 'creodias', 'onda', 'planetary_computer', 'sara']
Prostorový rozsah zájmového území musíme transformovat do WGS-84:
In [12]:
Copied!
extent = lucas_data.to_crs(epsg=4326).total_bounds
extent
extent = lucas_data.to_crs(epsg=4326).total_bounds
extent
Out[12]:
array([13.75920804, 50.41775 , 14.07773 , 50.55307 ])
V našem případě se budeme dotazovat služby Copernicus Data Space.
In [13]:
Copied!
dag.set_preferred_provider("cop_dataspace")
results = dag.search(
productType="S2_MSI_L2A",
start="2025-09-01",
end="2025-09-30",
geom=tuple(extent),
count=True
)
results.number_matched
dag.set_preferred_provider("cop_dataspace")
results = dag.search(
productType="S2_MSI_L2A",
start="2025-09-01",
end="2025-09-30",
geom=tuple(extent),
count=True
)
results.number_matched
2025-10-02 12:52:01,676 - eodag.core - INFO - core._do_search - Searching on provider cop_dataspace 2025-10-02 12:52:01,677 - eodag.search.base - INFO - base.get_sort_by_arg - cop_dataspace is configured with default sorting by 'startTimeFromAscendingNode' in ascending order 2025-10-02 12:52:01,685 - eodag.search.qssearch - INFO - qssearch._request - Sending search request: https://catalogue.dataspace.copernicus.eu/odata/v1/Products?$filter=Collection/Name%20eq%20%27SENTINEL-2%27%20and%20OData.CSC.Intersects%28area=geography%27SRID=4326%3BPOLYGON%20%28%2813.7592%2050.4177%2C%2013.7592%2050.5531%2C%2014.0777%2050.5531%2C%2014.0777%2050.4177%2C%2013.7592%2050.4177%29%29%27%29%20and%20Attributes/OData.CSC.StringAttribute/any%28att:att/Name%20eq%20%27productType%27%20and%20att/OData.CSC.StringAttribute/Value%20eq%20%27S2MSI2A%27%29%20and%20ContentDate/Start%20lt%202025-09-30T00:00:00.000Z%20and%20ContentDate/End%20gt%202025-09-01T00:00:00.000Z&$orderby=ContentDate/Start%20asc&$count=True&$top=20&$skip=0&$expand=Attributes&$expand=Assets 2025-10-02 12:52:38,958 - eodag.core - INFO - core._do_search - Found 16 result(s) on provider 'cop_dataspace'
Out[13]:
16
Snímky omezíme maximálním stupněm pokrytí oblačnosti.
In [14]:
Copied!
results_cc = results.filter_property(operator="le", cloudCover=10)
results_cc
results_cc = results.filter_property(operator="le", cloudCover=10)
results_cc
2025-10-02 12:52:38,963 - eodag.crunch.property - INFO - filter_property.proceed - Finished filtering products. 2 resulting products
Out[14]:
| SearchResult (2) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
0 EOProduct(id=S2A_MSIL2A_20250901T100701_N0511_R022_T33UVR_20250901T121415, provider=cop_dataspace)
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
1 EOProduct(id=S2A_MSIL2A_20250921T100731_N0511_R022_T33UVS_20250921T121814, provider=cop_dataspace)
|
Úkol: přidejte filtr, který vrátí pouze snímky, které kompletně pokrývají zájmové území.
Pro další kroky již budeme potřebovat nastavit přístupové údaje.
- Vytvoříme si učet: zde
- Vytvoříme konfigurační soubor
In [15]:
Copied!
username = "???"
password = "???"
dag.update_providers_config(f"""
cop_dataspace:
download:
extract: False
output_dir: /tmp/sentinel
delete_archive: False
auth:
credentials:
username: {username}
password: {password}
""")
username = "???"
password = "???"
dag.update_providers_config(f"""
cop_dataspace:
download:
extract: False
output_dir: /tmp/sentinel
delete_archive: False
auth:
credentials:
username: {username}
password: {password}
""")
Scény před stažením prohlédneme (zdroj):
In [16]:
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/S2A_MSIL2A_20250901T100701_N0511_R022_T33UVR_20250901T121415: 100%|█| 38.3k/38.3k [00:00<00:00, 339kB/s
2025-10-02 12:54:53,291 - eodag.product - INFO - _product.get_quicklook - Download recorded in /tmp/sentinel/quicklooks/S2A_MSIL2A_20250901T100701_N0511_R022_T33UVR_20250901T121415
quicklooks/S2A_MSIL2A_20250901T100701_N0511_R022_T33UVR_20250901T121415: 100%|█| 38.3k/38.3k [00:00<00:00, 334kB/s quicklooks/S2A_MSIL2A_20250921T100731_N0511_R022_T33UVS_20250921T121814: 0%| | 0.00/31.6k [00:00<?, ?B/s]
2025-10-02 12:54:53,932 - eodag.product - INFO - _product.get_quicklook - Download recorded in /tmp/sentinel/quicklooks/S2A_MSIL2A_20250921T100731_N0511_R022_T33UVS_20250921T121814
quicklooks/S2A_MSIL2A_20250921T100731_N0511_R022_T33UVS_20250921T121814: 100%|█| 31.6k/31.6k [00:00<00:00, 376kB/s
In [17]:
Copied!
results[1].properties['cloudCover']
results[1].properties['cloudCover']
Out[17]:
9.748491
Vybranou scénu nyní stáhneme:
In [19]:
Copied!
data_path = dag.download(results_cc[1])
data_path = dag.download(results_cc[1])
2025-10-02 17:07:39,276 - eodag.core - INFO - core.download - Local product detected. Download skipped
In [20]:
Copied!
from pathlib import Path
print(data_path, Path(data_path).exists())
from pathlib import Path
print(data_path, Path(data_path).exists())
/tmp/sentinel/S2A_MSIL2A_20250921T100731_N0511_R022_T33UVS_20250921T121814.SAFE.zip True
Přístup k datům¶
Vyzkoušíme si přístup k jednotlivým kanálům pomocí Rasterio.
In [21]:
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/sentinel/S2A_MSIL2A_20250921T100731_N0511_R022_T33UVS_20250921T121814.SAFE.zip/S2A_MSIL2A_20250921T100731_N0511_R022_T33UVS_20250921T121814.SAFE/MTD_MSIL2A.xml:10m:EPSG_32633
/home/martin/git/k155cvut/yusu/venv/lib/python3.13/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 [22]:
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
In [23]:
Copied!
with rio.open(sd[0]) as b10m:
b4 = b10m.read(1)
b4.shape
with rio.open(sd[0]) as b10m:
b4 = b10m.read(1)
b4.shape
Out[23]:
(10980, 10980)
In [24]:
Copied!
from rasterio.plot import show
show(b4, cmap=plt.cm.Greys_r, title="RED band")
from rasterio.plot import show
show(b4, cmap=plt.cm.Greys_r, title="RED band")
Out[24]:
<Axes: title={'center': 'RED band'}>
Ořez kanálů¶
Kanály ořežeme podle zájmového území.
In [25]:
Copied!
with rio.open(sd[0]) as b10m:
print(b10m.crs)
with rio.open(sd[0]) as b10m:
print(b10m.crs)
EPSG:32633
In [26]:
Copied!
from shapely.geometry import box
extent_utm = box(*lucas_data.to_crs(epsg=32633).total_bounds)
print(extent_utm)
from shapely.geometry import box
extent_utm = box(*lucas_data.to_crs(epsg=32633).total_bounds)
print(extent_utm)
POLYGON ((434627.3391108293 5585491.6076122215, 434627.3391108293 5600816.91057989, 411869.90084852505 5600816.91057989, 411869.90084852505 5585491.6076122215, 434627.3391108293 5585491.6076122215))
In [27]:
Copied!
from rasterio.mask import mask
with rio.open(sd[0]) as b10m:
clip_image, clip_transform = mask(b10m, [extent_utm], crop=True)
clip_meta = b10m.meta.copy()
clip_meta.update({
"height": clip_image.shape[1],
"width": clip_image.shape[2],
"transform": clip_transform
})
clip_image.shape
from rasterio.mask import mask
with rio.open(sd[0]) as b10m:
clip_image, clip_transform = mask(b10m, [extent_utm], crop=True)
clip_meta = b10m.meta.copy()
clip_meta.update({
"height": clip_image.shape[1],
"width": clip_image.shape[2],
"transform": clip_transform
})
clip_image.shape
Out[27]:
(6, 1062, 2277)
In [28]:
Copied!
show(clip_image[0, :, :], cmap=plt.cm.Greys_r, title="RED band")
show(clip_image[0, :, :], cmap=plt.cm.Greys_r, title="RED band")
Out[28]:
<Axes: title={'center': 'RED band'}>
Rasterizace trénovacích ploch¶
Trénovací plochu definujeme obalovou zónou 100m kolem LUCAS bodu.
In [29]:
Copied!
lucas_data.buffer(100)
lucas_data.buffer(100)
Out[29]:
0 POLYGON ((4594100.37 3042000.543, 4594099.888 ...
1 POLYGON ((4596100.851 3043999.936, 4596100.37 ...
2 POLYGON ((4600100.672 3048000.167, 4600100.191...
3 POLYGON ((4610100.128 3050000.214, 4610099.646...
4 POLYGON ((4610100.259 3053999.489, 4610099.778...
...
145 POLYGON ((4592100 3056000, 4592099.518 3055990...
146 POLYGON ((4608104.862 3050039.808, 4608104.381...
147 POLYGON ((4608099.183 3054001.125, 4608098.702...
148 POLYGON ((4606096.226 3050003.335, 4606095.744...
149 POLYGON ((4592104.383 3048025.627, 4592103.902...
Length: 150, dtype: geometry
Budeme uvažovat pouze první uroveň LUCAS nomenklatury:
In [30]:
Copied!
lucas_data["lc1_h"].apply(lambda x: x[0])
lucas_data["lc1_h"].apply(lambda x: x[0])
Out[30]:
0 B
1 B
2 B
3 B
4 B
..
145 C
146 F
147 D
148 A
149 A
Name: lc1_h, Length: 150, dtype: object
In [31]:
Copied!
# Písmeno 'A' odpovídá ASCII hodnotě 65
labels = lucas_data["lc1_h"].apply(lambda x: ord(x[0])-64)
labels
# Písmeno 'A' odpovídá ASCII hodnotě 65
labels = lucas_data["lc1_h"].apply(lambda x: ord(x[0])-64)
labels
Out[31]:
0 2
1 2
2 2
3 2
4 2
..
145 3
146 6
147 4
148 1
149 1
Name: lc1_h, Length: 150, dtype: int64
In [32]:
Copied!
from rasterio import features
clip_meta.update({
"driver": "GTiff",
"count": 1,
})
lucas_rast_fn = "/tmp/lucas_rast.tif"
with rio.open(lucas_rast_fn, "w+", **clip_meta) as out:
lucas_arr = out.read(1)
shapes = ((geom, value) for geom, value in zip(lucas_data.to_crs(epsg=32633).buffer(100), labels))
burned = features.rasterize(shapes=shapes, fill=0, out=lucas_arr, transform=out.transform)
out.write_band(1, burned)
from rasterio import features
clip_meta.update({
"driver": "GTiff",
"count": 1,
})
lucas_rast_fn = "/tmp/lucas_rast.tif"
with rio.open(lucas_rast_fn, "w+", **clip_meta) as out:
lucas_arr = out.read(1)
shapes = ((geom, value) for geom, value in zip(lucas_data.to_crs(epsg=32633).buffer(100), labels))
burned = features.rasterize(shapes=shapes, fill=0, out=lucas_arr, transform=out.transform)
out.write_band(1, burned)
In [33]:
Copied!
import numpy as np
np.unique(lucas_arr, return_counts=True)
import numpy as np
np.unique(lucas_arr, return_counts=True)
Out[33]:
(array([0, 1, 2, 3, 4, 5, 6, 7], dtype=uint16),
array([2405143, 1000, 5712, 2848, 944, 1894, 318,
315]))
In [34]:
Copied!
plt.imshow(lucas_arr)
plt.colorbar()
plt.imshow(lucas_arr)
plt.colorbar()
Out[34]:
<matplotlib.colorbar.Colorbar at 0x7fd22d5c27b0>
Operace překrytí¶
In [35]:
Copied!
plt.subplot(1, 3, 1)
plt.imshow(clip_image[0, :, :], cmap=plt.cm.Greys_r)
plt.title('Red band')
plt.subplot(1, 3, 2)
plt.imshow(clip_image[3, :, :], cmap=plt.cm.Greys)
plt.title('NIR band')
plt.subplot(1, 3, 3)
plt.imshow(lucas_arr, cmap=plt.cm.Spectral, interpolation='nearest')
plt.title('Training Data')
plt.show()
plt.subplot(1, 3, 1)
plt.imshow(clip_image[0, :, :], cmap=plt.cm.Greys_r)
plt.title('Red band')
plt.subplot(1, 3, 2)
plt.imshow(clip_image[3, :, :], cmap=plt.cm.Greys)
plt.title('NIR band')
plt.subplot(1, 3, 3)
plt.imshow(lucas_arr, cmap=plt.cm.Spectral, interpolation='nearest')
plt.title('Training Data')
plt.show()
Úkol: Zobrazte NIR ve zvýrazněných odstínech šedi podle histogramu.
In [67]:
Copied!
image_arr = np.dstack([clip_image[0, :, :], clip_image[3, :, :]])
image_arr.shape
image_arr.shape
np.unique(lucas_arr, return_counts=True)
image_arr = np.dstack([clip_image[0, :, :], clip_image[3, :, :]])
image_arr.shape
image_arr.shape
np.unique(lucas_arr, return_counts=True)
Out[67]:
(array([0, 1, 2, 3, 4, 5, 6, 7], dtype=uint16),
array([2405143, 1000, 5712, 2848, 944, 1894, 318,
315]))
In [69]:
Copied!
x = image_arr[lucas_arr > 0]
y = lucas_arr[lucas_arr > 0]
x.shape, y.shape # np.count_nonzero(lucas_arr > 0))
x = image_arr[lucas_arr > 0]
y = lucas_arr[lucas_arr > 0]
x.shape, y.shape # np.count_nonzero(lucas_arr > 0))
Out[69]:
((13031, 2), (13031,))
Úkol: Na závěr vykreslíme scatterplot, abychom zjistili, zda mezi proměnnými existuje lineární vztah.
In [39]:
Copied!
plt.scatter(x[:, 0], x[:, 1])
plt.scatter(x[:, 0], x[:, 1])
Out[39]:
<matplotlib.collections.PathCollection at 0x7fd22d5c30e0>
Závěrečné poznámky¶
- Pro zjednodušení používáme všechna LUCAS měření (z období 2006 až 2018) a Sentinel-2 z roku 2025. U reálného zadaní bychom stahovali satelitní snímky odpovídající časové značce LUCAS měření (atribut
survey_date). - Obalová zóna 100m je silně nadhodnocena z vizualizačních důvodů, nejde ale rozhodně o reprezentativní plochy. Obalová zóna by mohla odpovídat např. prostorovému rozlišení satelitních dat.