Dans ce tutoriel, nous ferons passer la structure de données par plusieurs formats, en lecture et en écriture : JSON, CSV, Parquet, GeoTIFF, GeoJSON, netCDF/HDF5 et Zarr. À la fin, nous comparons la taille sur disque des mêmes données selon les formats.
Tous les fichiers créés ou téléchargés dans ce carnet (notebook) vont dans un dossier local data/.
import os
os.makedirs('data', exist_ok=True)JSON (JavaScript Object Notation)¶
JSON est un format de données léger, lisible par un humain, utilisé dans les applications web et les API pour l’échange de données. On l’emploie pour stocker des métadonnées, des fichiers de configuration et de petits jeux de données, en particulier dans les applications web ou pour interroger des API (requêtes de données météo ou d’information géospatiale, par exemple). En voici un exemple simple :
{
"location": "Yellowstone",
"coordinates": {
"latitude": 44.423691,
"longitude": -110.588516
},
"elevation_m": 2399,
"temperature_c": 22.5
}
L’encodage des caractères est UTF-8. Les types de données d’un fichier JSON peuvent être : nombre, chaîne de caractères, booléen, tableau, objet (collection de paires nom-valeur) ou null. Plus d’informations sur JSON dans le cours EarthDataScience (en anglais).
En Python, on crée un fichier JSON avec la bibliothèque json :
import json
data = {
"location": "Yellowstone",
"coordinates": {"latitude": 44.423691, "longitude": -110.588516},
"elevation_m": 2399,
"temperature_c": 22.5
}
with open('data/data.json', 'w') as outfile:
json.dump(data, outfile)Formats de données tabulaires¶
Les données tabulaires sont courantes en géosciences. Le format CSV est un format commode, lisible par un humain, pour les petits jeux de données, qui stocke les données par lignes. Le format Parquet est lisible par la machine, organisé en colonnes et compressé, conçu pour les grands jeux de données.
CSV (Comma-Separated Values) :¶
CSV est un format simple et très répandu pour les données tabulaires, souvent utilisé en géosciences pour partager et stocker de petits jeux de données (échantillons de sol, relevés environnementaux…).
Il stocke les données en texte brut, ce qui le rend facile à lire mais inefficace pour les gros volumes. On le lit le plus souvent avec pandas.
import pandas as pd
# Build a small table and write it to a CSV file
df = pd.DataFrame({
'station': ['RAIN', 'HELN', 'BAKR'],
'temperature_c': [12.5, 10.1, 8.3],
'precipitation_mm': [3.0, 1.2, 0.0],
})
df.to_csv('data/output.csv', index=False)
# Read it back
df = pd.read_csv('data/output.csv')
dfParquet :¶
Parquet est un format de stockage binaire, en colonnes, optimisé pour l’efficacité, en particulier sur les grands jeux de données. Il est très utilisé dans les environnements big data (stockage d’imagerie satellitaire ou de sorties de modèles climatiques, par exemple).
# Writing the same table to a Parquet file
df.to_parquet('data/output.parquet', index=False)
# Reading the Parquet file
df_read = pd.read_parquet('data/output.parquet')
print(df_read) station temperature_c precipitation_mm
0 RAIN 12.5 3.0
1 HELN 10.1 1.2
2 BAKR 8.3 0.0
Données géospatiales¶
Les principaux formats de données géospatiales dans cette leçon :
- GeoTIFF : un standard de métadonnées qui embarque l’information de géoréférencement dans un fichier TIFF (Tagged Image File Format). Le Cloud-Optimized GeoTIFF (COG) en est une variante organisée pour l’accès distant.
- GeoJSON : un format pour encoder en JSON diverses structures de données géographiques (points, lignes, polygones).
import folium
import geopandas as gpd
import h5py
import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
import pooch
import rasterio
import xarray as xr
from folium.plugins import MarkerCluster
from rasterio.plot import show1. Données raster¶
1.1 rasterio pour lire un GeoTIFF¶
Une donnée raster est toute donnée pixélisée (ou maillée) où chaque pixel est associé à une position géographique précise. La valeur d’un pixel peut être continue (l’altitude, par exemple) ou catégorielle (l’occupation du sol).
Le paquet Python rasterio (documentation) lit des formats tels que GeoTIFF.
Voyez aussi les supports d’introduction d’EarthDataScience et les tutoriels du GeoHackweek.
Nous allons télécharger un raster de topographie en relief ombré depuis Natural Earth. Le fichier s’appelle HYP_50M_SR et arrive sous forme d’archive zip. Nous utilisons pooch pour le télécharger dans data/ et le décompresser en une seule étape.
files = pooch.retrieve(
url="https://naciscdn.org/naturalearth/50m/raster/HYP_50M_SR.zip",
known_hash="sha256:de6faaee29c8707764852a36c0973aea184e8ec94b6eebf21d061c923c8025de",
fname="HYP_50M_SR.zip",
path="./data",
processor=pooch.Unzip(extract_dir="HYP_50M_SR"),
)
tif_file = [f for f in files if f.endswith(".tif")][0]
print(tif_file)/Users/marinedenolle/Dropbox/CLASSES/ESS490/curriculum-book/book/Chapter2-DataManipulation/data/HYP_50M_SR/HYP_50M_SR.tif
Ouvrons maintenant le GeoTIFF avec rasterio.
elevation = rasterio.open(tif_file)Un jeu de données rasterio n’est pas un fichier netCDF : il n’a pas de .variables. Il stocke à la place une ou plusieurs bandes. Les attributs de base sont .count (nombre de bandes), .width et .height (dimensions du raster) et .crs (système de référence de coordonnées).
print("number of bands:", elevation.count)
print("width :", elevation.width)
print("height:", elevation.height)
print("crs :", elevation.crs)number of bands: 3
width : 10800
height: 5400
crs : EPSG:4326
elevation.indexes(1, 2, 3)Devinez-vous comment obtenir les types de données des bandes ?
# type below
elevation.dtypes('uint8', 'uint8', 'uint8')Et les limites du jeu de données :
elevation.boundsBoundingBox(left=-179.99999999999997, bottom=-89.99999999998201, right=179.99999999996405, top=90.0)print(elevation.transform * (0, 0)) # North West corner
print(elevation.transform * (elevation.width, elevation.height)) # South East corner(-179.99999999999997, 90.0)
(179.99999999996405, -89.99999999998201)
Comment interpréter les données : il y a trois bandes, une pour chacune des couleurs rouge, vert et bleu :
print(elevation.colorinterp[0])
print(elevation.colorinterp[1])
print(elevation.colorinterp[2])3
4
5
print(np.min(elevation.read(1)), np.max(elevation.read(1)))
print(np.min(elevation.read(2)), np.max(elevation.read(2)))
print(np.min(elevation.read(3)), np.max(elevation.read(3)))59 255
79 255
75 255
elevation.read(1) renvoie la première bande sous forme de tableau numpy 2D. Traçons-la seule, puis les trois bandes en image couleur.
band1 = elevation.read(1)
plt.imshow(band1, cmap='gray')
plt.colorbar(shrink=0.5, label='red band value')
plt.title('HYP_50M_SR, band 1')
plt.show()
image = elevation.read()
show(image)
<Axes: >1.2 Geopandas pour lire un GeoJSON¶
GeoJSON est un cas particulier de JSON qui stocke des objets géographiques (points, lignes, polygones) avec leurs attributs. Plutôt que d’en télécharger un, nous allons construire un petit GeoJSON de toutes pièces : une FeatureCollection de quatre volcans de la chaîne des Cascades. Cela montre explicitement la structure du format.
volcanoes_geojson = {
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"properties": {"name": "Mt Rainier", "elevation_m": 4392},
"geometry": {"type": "Point", "coordinates": [-121.7603, 46.8523]},
},
{
"type": "Feature",
"properties": {"name": "Mt St Helens", "elevation_m": 2549},
"geometry": {"type": "Point", "coordinates": [-122.1944, 46.1912]},
},
{
"type": "Feature",
"properties": {"name": "Mt Baker", "elevation_m": 3286},
"geometry": {"type": "Point", "coordinates": [-121.8144, 48.7768]},
},
{
"type": "Feature",
"properties": {"name": "Glacier Peak", "elevation_m": 3213},
"geometry": {"type": "Point", "coordinates": [-121.1132, 48.1125]},
},
],
}
with open('data/cascade_volcanoes.geojson', 'w') as f:
json.dump(volcanoes_geojson, f, indent=2)Notez la structure : chaque Feature a une geometry (ici un Point, longitude d’abord, latitude ensuite) et un dictionnaire properties pour les attributs. Relisez maintenant le fichier avec geopandas :
volcanoes = gpd.read_file('data/cascade_volcanoes.geojson')
volcanoesGeopandas a interprété la colonne de géométrie et fixé le système de référence de coordonnées à WGS84 (EPSG:4326), le défaut de GeoJSON :
volcanoes.crs<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: GreenwichTraçons les données sur une carte interactive. Folium est un agréable paquet Python de visualisation. Le tutoriel Folium du Geohackweek est lui aussi instructif.
m = folium.Map(location=[47.5, -121.8], zoom_start=7)
marker_cluster = MarkerCluster().add_to(m)
for _, row in volcanoes.iterrows():
folium.Marker(
location=[row.geometry.y, row.geometry.x],
popup=f"{row['name']}: {row['elevation_m']} m",
).add_to(marker_cluster)
m2. Formats hiérarchiques : NETCDF4 & HDF5¶
Les formats de données hiérarchiques sont conçus pour stocker de grands volumes de données dans un seul fichier. Ils imitent un système de fichiers (une arborescence de répertoires imbriqués) à l’intérieur d’un fichier unique. Deux formats hiérarchiques dominent (HDF5 et NETCDF4), plus un conçu pour le cloud (Zarr). Les formats hiérarchiques peuvent stocker de nombreux types de données (numériques vs chaînes de caractères).
HDF5¶
Le Hierarchical Data Format version 5 (HDF5) est un format de fichier open source qui supporte des données volumineuses, complexes et hétérogènes. HDF5 utilise une structure de type « répertoire de fichiers » qui permet d’organiser les données à l’intérieur du fichier de multiples façons structurées, comme vous le feriez avec les fichiers de votre ordinateur. Le format HDF5 permet aussi d’embarquer des métadonnées, ce qui le rend auto-descriptif. Tous les éléments (le fichier lui-même, les groupes et les jeux de données) peuvent porter des métadonnées qui décrivent l’information contenue dans l’élément.
Exemple de structure HDF :
- les jeux de données (datasets), qui sont des tableaux multidimensionnels typés ;
- les groupes, structures conteneurs qui peuvent contenir des jeux de données et d’autres groupes.
Figure : exemple de données HDF5, tiré de [neonscience](https://www.neonscience.org/resources/learning-hub/tutorials/about-hdf5)NetCDF¶
Le network Common Data Form, ou netCDF, a été créé au début des années 1990 pour résoudre certaines des difficultés du travail avec des tableaux N-dimensionnels. NetCDF est un ensemble de formats de données binaires auto-descriptifs, indépendants de la machine, et d’outils logiciels qui facilitent la création, l’accès et le partage de données scientifiques stockées en tableaux N-dimensionnels, accompagnées des métadonnées décrivant le contenu de chaque tableau. NetCDF a été bâti par la communauté des sciences du climat à une époque où les modèles climatiques régionaux commençaient à produire des fichiers de sortie de plus en plus gros. NetCDF version 4 est un sous-ensemble de HDF5 : les fichiers netCDF4 s’ouvrent donc avec les outils HDF5.
Manipuler de grands tableaux¶
Les formats netCDF et HDF5 n’imposent aucune limite de taille de fichier. En revanche, tout outil d’analyse qui charge en mémoire les données d’un tableau netCDF pour une opération de calcul est limité par la mémoire disponible de la machine.
Mais lents en entrées-sorties¶
À la lecture d’un fichier hiérarchique, tout l’arbre de la structure de données est parcouru depuis la racine. Comme ce parcours recommence à chaque requête de l’utilisateur, la lecture de HDF5 et de netCDF est lente comparée aux formats conçus pour l’accès parallèle ou dans le cloud.
2.1 Créer un jeu de données et l’écrire en netCDF¶
Nous créons un champ synthétique d’anomalies mensuelles de température avec le paquet du cours mlgeo_synth, gardons une tranche temporelle comme champ 2D, et l’enveloppons dans un Dataset xarray avec coordonnées et unités.
import mlgeo_synth
field, truth = mlgeo_synth.climate_field(n_lat=40, n_lon=80, n_months=120, seed=42)
print(field.shape) # (n_months, n_lat, n_lon)(120, 40, 80)
ds = xr.Dataset(
data_vars={
"temperature_anomaly": (
("lat", "lon"),
field[0],
{"units": "degC", "long_name": "monthly temperature anomaly"},
)
},
coords={
"lat": ("lat", truth["lat"], {"units": "degrees_north"}),
"lon": ("lon", truth["lon"], {"units": "degrees_east"}),
},
attrs={"title": "Synthetic temperature anomaly, month 0", "source": "mlgeo_synth.climate_field"},
)
dsds.to_netcdf('data/temperature_anomaly.nc')2.2 Relire le fichier avec xarray¶
ds_read = xr.open_dataset('data/temperature_anomaly.nc')
ds_read.temperature_anomaly.plot()
plt.title('Synthetic temperature anomaly (month 0)')
plt.show()
ds_read.close()
2.3 Le même fichier à travers les bibliothèques netCDF4 et h5py¶
La bibliothèque netCDF4 montre la structure auto-descriptive : dimensions, variables et attributs.
geo = nc.Dataset('data/temperature_anomaly.nc')
print(geo)
print(geo['temperature_anomaly'])
geo.close()<class 'netCDF4.Dataset'>
root group (NETCDF4 data model, file format HDF5):
title: Synthetic temperature anomaly, month 0
source: mlgeo_synth.climate_field
dimensions(sizes): lat(40), lon(80)
variables(dimensions): float64 temperature_anomaly(lat, lon), float64 lat(lat), float64 lon(lon)
groups:
<class 'netCDF4.Variable'>
float64 temperature_anomaly(lat, lon)
_FillValue: nan
units: degC
long_name: monthly temperature anomaly
unlimited dimensions:
current shape = (40, 80)
filling on
Comme les fichiers netCDF4 sont des fichiers HDF5, h5py peut ouvrir le même fichier. Les variables apparaissent comme des jeux de données HDF5, et les métadonnées comme des attributs HDF5.
with h5py.File('data/temperature_anomaly.nc', 'r') as f:
print("datasets in the root group:", list(f.keys()))
dset = f['temperature_anomaly']
print("shape:", dset.shape, "dtype:", dset.dtype)
print("attributes:", dict(dset.attrs))datasets in the root group: ['temperature_anomaly', 'lat', 'lon']
shape: (40, 80) dtype: float64
attributes: {'_Netcdf4Coordinates': array([0, 1], dtype=int32), '_FillValue': array([nan]), 'units': np.bytes_(b'degC'), 'long_name': np.bytes_(b'monthly temperature anomaly'), 'DIMENSION_LIST': array([array([<HDF5 object reference>], dtype=object),
array([<HDF5 object reference>], dtype=object)], dtype=object), '_Netcdf4Dimid': np.int32(0)}
3. Zarr¶
Zarr est un format optimisé pour le cloud, destiné aux tableaux N-dimensionnels. Au lieu d’un fichier unique, un store Zarr est un répertoire de blocs (chunks) compressés, accompagnés de petits fichiers de métadonnées JSON. Les stockages objet du cloud et les tâches parallèles peuvent lire des blocs individuels sans parcourir toute l’arborescence, ce qui supprime le goulot d’étranglement d’entrées-sorties de netCDF/HDF5.
Xarray écrit directement en Zarr :
ds.to_zarr('data/temperature_anomaly.zarr', mode='w')/Users/marinedenolle/Dropbox/CLASSES/ESS490/curriculum-book/.pixi/envs/default/lib/python3.12/site-packages/zarr/api/asynchronous.py:246: ZarrUserWarning: Consolidated metadata is currently not part in the Zarr format 3 specification. It may not be supported by other zarr implementations and may change in the future.
warnings.warn(
<xarray.backends.zarr.ZarrStore at 0x38823f4c0>ds_zarr = xr.open_zarr('data/temperature_anomaly.zarr')
ds_zarrL’aller-retour préserve les données et les métadonnées :
np.allclose(ds_zarr.temperature_anomaly.values, ds.temperature_anomaly.values)True4. Comparer les tailles sur disque selon les formats¶
Nous stockons maintenant les mêmes données dans les formats ci-dessus et comparons leurs tailles sur disque.
- Pour la comparaison table, nous aplatissons le champ 2D en une table de données (data frame, une ligne par cellule de grille) et l’écrivons en CSV et en Parquet.
- Pour la comparaison tableau, nous réutilisons le fichier netCDF et le store Zarr déjà écrits. Un store Zarr étant un répertoire, nous additionnons les tailles de tous les fichiers qu’il contient.
# Flatten the array to a table: one row per (lat, lon) grid cell
table = ds.temperature_anomaly.to_dataframe().reset_index()
table.to_csv('data/temperature_anomaly.csv', index=False)
table.to_parquet('data/temperature_anomaly.parquet', index=False)
table.head()def path_size_bytes(path):
"""Size of a file, or the total size of all files under a directory."""
if os.path.isfile(path):
return os.path.getsize(path)
total = 0
for root, _, filenames in os.walk(path):
for filename in filenames:
total += os.path.getsize(os.path.join(root, filename))
return total
sizes = {
'CSV (table)': path_size_bytes('data/temperature_anomaly.csv'),
'Parquet (table)': path_size_bytes('data/temperature_anomaly.parquet'),
'netCDF (array)': path_size_bytes('data/temperature_anomaly.nc'),
'Zarr (array)': path_size_bytes('data/temperature_anomaly.zarr'),
}
comparison = pd.DataFrame(
{'size_kB': [v / 1024 for v in sizes.values()]},
index=sizes.keys(),
).round(1)
print(comparison) size_kB
CSV (table) 133.0
Parquet (table) 33.0
netCDF (array) 33.9
Zarr (array) 29.9
Interprétation : le fichier CSV est le plus gros parce qu’il stocke chaque nombre sous forme de texte, chiffre par chiffre. Parquet stocke la même table dans une représentation binaire compressée et est plusieurs fois plus petit. Le fichier netCDF et le store Zarr contiennent le même tableau sous forme binaire et se retrouvent dans la même gamme de tailles ; à cette échelle, l’écart entre les deux vient du surcoût du conteneur (en-têtes HDF5 pour netCDF, fichiers de blocs et métadonnées JSON pour Zarr), pas des données elles-mêmes. L’écart entre formats texte et binaires grandit avec la taille du jeu de données, et la structure en blocs de Zarr devient payante quand de nombreux processus lisent en même temps différents blocs d’un même grand tableau.