Masquage de la qualité des nuages et des pixels pour Landsat
Produits utilisés : ls8_sr
Mots clés méthodes de données ; dilater, méthodes de données ; masquage, masquage des nuages, méthodes de données ; mise en mémoire tampon
Aperçu
Par le passé, les chercheurs en télédétection rejetaient les scènes partiellement couvertes de nuages au profit de scènes sans nuages. Cependant, les techniques d’analyse multitemporelle utilisent de plus en plus chaque pixel de qualité garantie dans une série temporelle d’observations. La capacité d’exclure automatiquement les pixels de mauvaise qualité (par exemple les nuages, les ombres ou d’autres données non valides) est essentielle pour permettre ces nouvelles méthodes de télédétection.
Les données satellite prêtes à être analysées de Digital Earth Africa incluent des informations sur la qualité des pixels qui peuvent être utilisées pour « masquer » facilement les données (c’est-à-dire conserver uniquement certains pixels dans une image) afin d’obtenir une série chronologique contenant uniquement des pixels clairs ou sans nuages.
Description
Dans ce bloc-notes, nous montrons comment masquer les données du satellite Digital Earth Africa Landsat à l’aide de masques booléens. Le bloc-notes montre comment :
Charger une série chronologique de données satellite incluant la bande de qualité de pixel « pixel_qa »
Inspectez les attributs « flags_definition » du groupe
Créez un masque où les pixels sont nuageux, ont une ombre de nuage ou ne contiennent aucune donnée
Appliquer des opérateurs morphologiques binaires sur les pixels nuageux pour améliorer le masque
Masquage des aérosols à forte concentration
Appliquer les masques aux données satellites afin de ne conserver que les observations de bonne qualité et tracer les résultats
Redimensionnement des données Landsat et masquage des données invalides
Utilisation de « load_ard » pour masquer les pixels de mauvaise qualité, en profitant de toutes ses fonctionnalités intégrées pour le masquage
Digital Earth Africa fournit des fonctions wrapper qui fourniront automatiquement des données satellite masquées par les nuages. Plus d’informations peuvent être trouvées dans le bloc-notes Using_load_ard.
Commencer
Nous importons d’abord les packages pertinents et nous nous connectons au datacube. Ensuite, nous définissons notre exemple de zone d’intérêt et chargeons une série chronologique de données satellite.
[1]:
%matplotlib inline
import datacube
from pprint import pprint
import geopandas as gpd
from datacube.utils.geometry import Geometry
from datacube.utils import masking
from odc.algo import mask_cleanup, erase_bad, to_f32
from deafrica_tools.plotting import rgb
from deafrica_tools.datahandling import mostcommon_crs, load_ard
from deafrica_tools.areaofinterest import define_area
Se connecter au datacube
[2]:
dc = datacube.Datacube(app="Masking_data")
Créer une requête et charger des données satellite
Pour démontrer comment masquer les données satellite, nous allons charger les données RVB de réflectance de surface Landsat 8 avec une bande de classification de la qualité des pixels appelée « pixel_quality ».
Pour définir la zone d’intérêt, deux méthodes sont disponibles :
En spécifiant la latitude, la longitude et la zone tampon. Cette méthode nécessite que vous saisissiez la latitude centrale, la longitude centrale et la valeur de la zone tampon en degrés carrés autour du point central que vous souhaitez analyser. Par exemple, « lat = 10,338 », « lon = -1,055 » et « buffer = 0,1 » sélectionneront une zone avec un rayon de 0,1 degré carré autour du point avec les coordonnées (10,338, -1,055).
By uploading a polygon as a
GeoJSON or Esri Shapefile. If you choose this option, you will need to upload the geojson or ESRI shapefile into the Sandbox using Upload Files button
in the top left corner of the Jupyter Notebook interface. ESRI shapefiles must be uploaded with all the related files (.cpg, .dbf, .shp, .shx). Once uploaded, you can use the shapefile or geojson to define the area of interest. Remember to update the code to call the file you have uploaded.
Pour utiliser l’une de ces méthodes, vous pouvez décommenter la ligne de code concernée et commenter l’autre. Pour commenter une ligne, ajoutez le symbole "#" avant le code que vous souhaitez commenter. Par défaut, la première option qui définit l’emplacement à l’aide de la latitude, de la longitude et du tampon est utilisée.
Si vous exécutez le bloc-notes pour la première fois, conservez les paramètres par défaut ci-dessous. Cela permettra de démontrer le fonctionnement de l’analyse et de fournir des résultats significatifs.
[3]:
# Method 1: Specify the latitude, longitude, and buffer
aoi = define_area(lat=35.7718, lon=-5.8096, buffer=0.03)
# Method 2: Use a polygon as a GeoJSON or Esri Shapefile.
# aoi = define_area(vector_path='aoi.shp')
#Create a geopolygon and geodataframe of the area of interest
geopolygon = Geometry(aoi["features"][0]["geometry"], crs="epsg:4326")
geopolygon_gdf = gpd.GeoDataFrame(geometry=[geopolygon], crs=geopolygon.crs)
# Get the latitude and longitude range of the geopolygon
lat_range = (geopolygon_gdf.total_bounds[1], geopolygon_gdf.total_bounds[3])
lon_range = (geopolygon_gdf.total_bounds[0], geopolygon_gdf.total_bounds[2])
# Create a reusable query
query = {
'x': lon_range,
'y': lat_range,
'time': ('2016-09', '2016-10'),
}
# Identify the most common projection system in the input query
output_crs = mostcommon_crs(dc=dc, product='ls8_sr', query=query)
# Load data from the Landsat-8
data = dc.load(product="ls8_sr",
measurements=["blue", "green", "red",
"pixel_quality", "qa_aerosol"],
output_crs=output_crs,
resolution=(-30, 30),
align=(15, 15),
**query)
print(data)
/home/jovyan/dev/deafrica-sandbox-notebooks/Tools/deafrica_tools/datahandling.py:733: UserWarning: Multiple UTM zones ['epsg:32629', 'epsg:32630'] were returned for this query. Defaulting to the most common zone: epsg:32630
warnings.warn(
<xarray.Dataset>
Dimensions: (time: 7, y: 228, x: 188)
Coordinates:
* time (time) datetime64[ns] 2016-09-02T11:03:08.816501 ... 2016-...
* y (y) float64 3.966e+06 3.966e+06 ... 3.959e+06 3.959e+06
* x (x) float64 2.432e+05 2.432e+05 ... 2.488e+05 2.488e+05
spatial_ref int32 32630
Data variables:
blue (time, y, x) uint16 8398 8338 8295 8216 ... 16903 19403 21435
green (time, y, x) uint16 8813 8771 8707 8622 ... 20944 23082 23060
red (time, y, x) uint16 8235 8199 8136 8077 ... 22769 24061 24180
pixel_quality (time, y, x) uint16 23888 23888 23888 ... 21824 21824 21824
qa_aerosol (time, y, x) uint8 224 224 224 224 224 ... 194 224 224 64 96
Attributes:
crs: epsg:32630
grid_mapping: spatial_ref
L’absence d’observation par satellite est indiquée par une valeur « nodata » pour la bande, qui est répertoriée dans la catégorie Attributs du « xarray.DataArray » renvoyé.
[4]:
data.red.attrs
[4]:
{'units': '1', 'nodata': 0, 'crs': 'epsg:32630', 'grid_mapping': 'spatial_ref'}
Nous voyons que l’attribut « nodata » renvoie la valeur « 0 ».
Nous pouvons retrouver le schéma de classification de la bande « pixel_qa » dans sa définition de drapeaux.
[5]:
flags_def = data.pixel_quality.attrs["flags_definition"]
pprint(flags_def)
{'cirrus': {'bits': 2,
'values': {'0': 'not_high_confidence', '1': 'high_confidence'}},
'cirrus_confidence': {'bits': [14, 15],
'values': {'0': 'none',
'1': 'low',
'2': 'reserved',
'3': 'high'}},
'clear': {'bits': 6, 'values': {'0': False, '1': True}},
'cloud': {'bits': 3,
'values': {'0': 'not_high_confidence', '1': 'high_confidence'}},
'cloud_confidence': {'bits': [8, 9],
'values': {'0': 'none',
'1': 'low',
'2': 'medium',
'3': 'high'}},
'cloud_shadow': {'bits': 4,
'values': {'0': 'not_high_confidence',
'1': 'high_confidence'}},
'cloud_shadow_confidence': {'bits': [10, 11],
'values': {'0': 'none',
'1': 'low',
'2': 'reserved',
'3': 'high'}},
'dilated_cloud': {'bits': 1, 'values': {'0': 'not_dilated', '1': 'dilated'}},
'nodata': {'bits': 0, 'values': {'0': False, '1': True}},
'snow': {'bits': 5,
'values': {'0': 'not_high_confidence', '1': 'high_confidence'}},
'snow_ice_confidence': {'bits': [12, 13],
'values': {'0': 'none',
'1': 'low',
'2': 'reserved',
'3': 'high'}},
'water': {'bits': 7, 'values': {'0': 'land_or_cloud', '1': 'water'}}}
Nous voyons que « pixel_quality » signale également les pixels « nodata » et, avec les pixels « cloud » et « cloud_shadow », il capte également les pixels « snow » et « water ».
Création d’un masque de qualité nuage et pixel
Nous créons un masque en spécifiant les conditions que nos pixels doivent satisfaire. Mais nous n’aurons besoin que des étiquettes (et non des valeurs) pour créer un masque.
Dans cet exemple ci-dessous, nous souhaitons exclure plusieurs catégories de pixels défectueux, c’est-à-dire que si un pixel a l’un des indicateurs (cloud, cirrus ou cloud_shadow) défini, le pixel sera masqué.
[6]:
quality_flags = dict(
cloud="high_confidence", # True where there is cloud
cirrus="high_confidence",# True where there is cirrus cloud
cloud_shadow="high_confidence",# True where there is cloud shadow
)
# set cloud_mask: True=cloud, False=non-cloud
mask, _= masking.create_mask_value(flags_def, **quality_flags)
#add the cloud mask to our dataset
data['cloud_mask'] = (data['pixel_quality'] & mask) != 0
Ci-dessous, nous allons tracer le masque de qualité des pixels avec les images satellites en couleurs réelles.
Le masque de nuages correspond-il exactement aux nuages que vous voyez dans les tracés RVB ? L’algorithme de qualité des pixels de Landsat présente des limites connues qui conduisent à classer par erreur des objets lumineux, tels que des plages et des villes, comme des nuages.
[7]:
#plot the locations where there are clouds and cloud shadows
data['cloud_mask'].plot(col="time", col_wrap=7, add_colorbar=False);
[8]:
# Plot in True colour
rgb(data, col="time", col_wrap=7)
Application du traitement morphologique sur le masque de nuages
Nous pouvons améliorer les faux positifs détectés par le masque de qualité des pixels de Landsat en appliquant des techniques de traitement d’image morphologique binaire (par exemple, binary_closing, binary_erosion etc.). La bibliothèque Open-Data-Cube odc-algo possède une fonction, odc.algo.mask_cleanup qui peut effectuer quelques-unes de ces opérations. Ci-dessous, nous allons essayer d’améliorer le masque de nuage en appliquant un certain nombre de ces filtres.
N’hésitez pas à expérimenter avec les valeurs dans « filtres »
[9]:
# set the filters to apply. The integers refer to the number of pixels
filters = [("opening", 2),("dilation", 2)]
[10]:
# Use the mask_cleanup function to apply the filters
data['cloud_mask_filtered'] = mask_cleanup(data['cloud_mask'], mask_filters=filters)
[11]:
#plot the results
data['cloud_mask_filtered'].plot(col="time", col_wrap=7, add_colorbar=False);
Masquage des pixels avec un aérosol à haute teneur
L’algorithme de correction atmosphérique appliqué à Landsat utilise l’aérosol comme entrée. Lorsque la valeur de l’aérosol est élevée, la valeur de la réflectance de surface peut être peu fiable, en particulier sur les surfaces sombres. Par conséquent, pour certaines applications, le masquage des données avec une valeur élevée de l’aérosol peut aider à éliminer les observations peu fiables.
Dans l’exemple ci-dessous, nous utilisons la fonction « masking.make_mask() » pour créer un masque pour les pixels avec « aérosol élevé » et le combiner avec le masque de nuage filtré ci-dessus.
[12]:
# find the bit mask flags and print them
aerosol_flags = data.qa_aerosol.attrs["flags_definition"]
pprint(aerosol_flags)
#create a dict of flags we want to mask
aerosol_quality_flags = dict(
aerosol_level="high",
)
#create a boolean mask where high-aerosol = True
data['aerosol_mask'] = masking.make_mask(data['qa_aerosol'], **aerosol_quality_flags)
{'aerosol_level': {'bits': [6, 7],
'values': {'0': 'climatology',
'1': 'low',
'2': 'medium',
'3': 'high'}},
'interpolated_aerosol': {'bits': 5, 'values': {'0': False, '1': True}},
'nodata': {'bits': 0, 'values': {'0': False, '1': True}},
'valid_aerosol_retrieval': {'bits': 1, 'values': {'0': False, '1': True}},
'water': {'bits': 2, 'values': {'0': False, '1': True}}}
[13]:
data['aerosol_mask'].plot(col="time", col_wrap=7);
Nous combinons maintenant le masque aérosol avec le masque « clear_filtered » pour produire un masque où les données ne sont pas fiables en raison de nuages ou d’aérosols élevés
[14]:
data['cloud_aerosol_mask'] = data['aerosol_mask'] | data['cloud_mask_filtered']
#plot the result
data['cloud_aerosol_mask'].plot(col="time", col_wrap=7);
Application du masque nuage
Nous pouvons maintenant obtenir les images claires que nous souhaitons en effaçant les nuages, les pixels sans données et/ou les pixels à forte concentration d’aérosols des données. Nous procéderons ainsi pour le masque PQ d’origine, le masque PQ qui a subi un traitement d’image morphologique binaire et le masque avec masquage supplémentaire des aérosols.
[15]:
# erase pixels with cloud
clear = erase_bad(data.drop_vars(['cloud_mask_filtered', 'cloud_mask', 'cloud_aerosol_mask', 'aerosol_mask']),
data['cloud_mask'])
# erase pixels with the cloud_filtering
clear_filtered = erase_bad(data.drop_vars(['cloud_mask_filtered', 'cloud_mask', 'cloud_aerosol_mask','aerosol_mask']),
data['cloud_mask_filtered'])
# erase pixels with cloud and high aerosol
clear_filtered_aerosol = erase_bad(data.drop_vars(['cloud_mask_filtered', 'cloud_mask', 'cloud_aerosol_mask','aerosol_mask']),
data['cloud_aerosol_mask'])
Tracez les résultats de notre masquage « clair »
[16]:
rgb(clear, col="time", col_wrap=7)
Tracez les résultats de notre masquage « clear_filtered »
Comme vous pouvez le constater, les opérations de filtrage morphologique ont minimisé l’impact des faux positifs dans le masque nuageux au-dessus des villes et le long de la côte
[17]:
rgb(clear_filtered, col="time", col_wrap=7)
Tracez les résultats de nos données masquées « clear_filtered » et « aerosol »
Le masquage des aérosols à forte concentration supprime de nombreux pixels au-dessus de l’océan où les valeurs de réflectance de surface sont faibles et peu fiables. Il faut cependant noter que cela peut également supprimer des observations utiles. Les utilisateurs doivent utiliser ce masque avec précaution.
[18]:
rgb(clear_filtered_aerosol, col="time", col_wrap=7)
Rappel des données Landsat et masquage des données invalides
Redimensionnez les données Landsat à ~0-1 conformément aux facteurs d’échelle et de décalage fournis par la documentation de la collection USGS Landsat 2 <https://docs.digitalearthafrica.org/en/latest/data_specs/Landsat_C2_SR_specs.html#Measurements>`__.
Remarque : les données USGS Landsat présentent des valeurs d’échelle et de décalage différentes selon la mesure. Les exemples d’échelle et de valeurs de décalage indiqués ci-dessous s’appliquent uniquement aux bandes de réflectance de surface
Avant de pouvoir redimensionner les valeurs, nous devons nous assurer que les données sont d’abord converties en « float32 », ce qui garantit que les valeurs sans données dans le type de données natif (« 0 ») sont converties en « NaN », garantissant ainsi que les valeurs sans données dans les données redimensionnées sont également « NaN ».
[19]:
# Convert the data to `float32` and set no-data values to `NaN`:
clear_filtered = to_f32(clear_filtered)
[20]:
#Apply the scale and offset factors to the Landsat data
sr_bands = ["blue", "green", "red"]
for band in sr_bands:
clear_filtered[band] = 2.75e-5 * clear_filtered[band] - 0.2
Ci-dessous, nous retraçons « clear_filtered » en vraies couleurs. Notez que la surface du terrain semble plus vive car la fonction « rgb » n’inclut plus les valeurs « 0 » dans les rampes de couleurs et parce que la réflectance de surface a été rééchelonnée à ses valeurs prévues (~ 0-1).
[21]:
rgb(clear_filtered, col="time", col_wrap=7)
/usr/local/lib/python3.10/dist-packages/matplotlib/cm.py:478: RuntimeWarning: invalid value encountered in cast
xx = (xx * 255).astype(np.uint8)
Masquage des nuages, filtrage morphologique et redimensionnement avec load_ard
La plupart des étapes décrites ci-dessus peuvent être réalisées à l’aide de la fonction « deafrica_tools.datahanding.load_ard ». Cette fonction créera un masque de nuages, appliquera des filtres morphologiques (si vous le souhaitez) et appliquera les facteurs d’échelle et de décalage Landsat.
Remarque : « load_ard » présente certaines limitations en raison de sa conception. Une partie importante de sa fonctionnalité prévue est de fournir un moyen simple de concaténer Landsat 5, 7 et 8 ensemble pour former un seul objet xarray. Cependant, cela signifie que « load_ard » ne peut appliquer un masquage de qualité de pixel qu’aux catégories pq qui sont communes à tous les capteurs Landsat. Par exemple, Landsat-8 dispose d’un indicateur de bits dédié pour les bandes de cirrus, mais Landsat 5 et 7 ne l’ont pas. Cela signifie que « load_ard » ne peut pas accepter « “cirrus”: “high_confidence” » comme catégorie pq. Le même problème s’applique également au masquage des aérosols, Landat 5 et 7 ont des moyens différents d’enregistrer les aérosols à forte concentration que Landsat 8, et donc « load_ard » ne prend pas en charge le masquage des aérosols. Cependant, « load_ard » supprime les valeurs de réflectance de surface négatives (les valeurs négatives sont un artefact occasionnel des algorithmes de correction employés par l’USGS), qui coïncident souvent avec des zones sujettes à des valeurs d’aérosols élevées telles que les grands plans d’eau et le long de la côte.
[22]:
# set the filters to apply
filters = [("opening", 2),("dilation", 2)]
# # Load data from the Landsat-8
data = load_ard(dc=dc,
products=["ls8_sr"],
measurements=["blue", "green", "red"],
mask_filters=filters,
output_crs=output_crs,
resolution=(-30, 30),
align=(15,15),
**query)
# plot the data
rgb(data, col="time", col_wrap=7)
Using pixel quality parameters for USGS Collection 2
Finding datasets
ls8_sr
Applying morphological filters to pq mask [('opening', 2), ('dilation', 2)]
Applying pixel quality/cloud mask
Re-scaling Landsat C2 data
Loading 7 time steps
/usr/local/lib/python3.10/dist-packages/matplotlib/cm.py:478: RuntimeWarning: invalid value encountered in cast
xx = (xx * 255).astype(np.uint8)
Informations Complémentaires
Licence : Le code de ce carnet est sous licence Apache, version 2.0 <https://www.apache.org/licenses/LICENSE-2.0>. Les données de Digital Earth Africa sont sous licence Creative Commons par attribution 4.0 <https://creativecommons.org/licenses/by/4.0/>.
Contact : Si vous avez besoin d’aide, veuillez poster une question sur le canal Slack Open Data Cube <http://slack.opendatacube.org/>`__ ou sur le GIS Stack Exchange en utilisant la balise open-data-cube (vous pouvez consulter les questions posées précédemment ici). Si vous souhaitez signaler un problème avec ce bloc-notes, vous pouvez en déposer un sur Github.
Version de Datacube compatible :
[23]:
print(datacube.__version__)
1.8.15
Dernier test :
[24]:
from datetime import datetime
datetime.today().strftime('%Y-%m-%d')
[24]:
'2023-08-11'