Rastériser des vecteurs et vectoriser des rasters

Mots-clés données utilisées ; WOfS, méthodes de données ; rastériser, méthodes de données ; vectoriser, format de données ; GeoTIFF, format de données ; shapefile

Aperçu

De nombreux flux de travail de télédétection et/ou de géospatiale nécessitent une conversion entre des données vectorielles (par exemple, des fichiers de formes) et des données raster (par exemple, des données basées sur des pixels comme celles d’un « xarray.DataArray »). Par exemple, nous pouvons avoir besoin d’utiliser un fichier de formes comme masque pour limiter l’étendue de l’analyse d’un raster, ou avoir des données raster que nous souhaitons convertir en données vectorielles pour permettre des opérations géométriques faciles.

Description

Dans ce bloc-notes, nous montrons comment utiliser la fonction Digital Earth Africa xr_rasterize et xr_vectorize dans deafrica_tools.spatial. Le bloc-notes montre comment :

  1. Charger les données du produit Water Observations from Space (WOfS)

  2. Vectorisez l’objet WOfS « xarray.DataArray » basé sur des pixels en un objet « geopandas.GeoDataFrame » basé sur des vecteurs contenant des plans d’eau persistants sous forme de polygones

  3. Exporter le « geopandas.GeoDataFrame » en tant que fichier de formes

  4. Rastérisez les données vectorielles « geopandas.GeoDataFrame » dans un objet « xarray.DataArray » et exportez les résultats au format GeoTIFF

Commencer

Pour exécuter cette analyse, exécutez toutes les cellules du bloc-notes, en commençant par la cellule « Charger les packages ».

Charger des paquets

[1]:

%matplotlib inline

import datacube
import geopandas as gpd
from odc.geo.geom import Geometry

from deafrica_tools.datahandling import mostcommon_crs
from deafrica_tools.spatial import xr_vectorize, xr_rasterize
from deafrica_tools.areaofinterest import define_area

Se connecter au datacube

[2]:

dc = datacube.Datacube(app='Rasterise_vectorise')

Charger les données WOfS à partir du datacube

Nous allons charger un résumé annuel du produit Observations de l’eau depuis l’espace (WOfS) pour nous fournir des données avec lesquelles travailler.

Pour définir la zone d’intérêt, deux méthodes sont disponibles :

  1. 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).

  2. 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 5d2be07e35344b0a8b4ba5a6c2d2db7e 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.

[3]:

# Select a location
# Method 1: Specify the latitude, longitude, and buffer
aoi = define_area(lat=13.50, lon=-15.42, buffer=0.2)

# 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': ('2017')
}
# Identify the most common projection system in the input query
output_crs = mostcommon_crs(dc=dc, product='ls8_sr', query=query)

# Load WoFS through the datacube
ds = dc.load(product='wofs_ls_summary_annual',
             output_crs=output_crs,
             align=(15, 15),
             resolution=(-30, 30),
             **query)

print(ds)

<xarray.Dataset> Size: 17MB
Dimensions:      (time: 1, y: 1478, x: 1447)
Coordinates:
  * time         (time) datetime64[ns] 8B 2017-07-02T11:59:59.999999
  * y            (y) float64 12kB 1.515e+06 1.515e+06 ... 1.47e+06 1.47e+06
  * x            (x) float64 12kB 4.328e+05 4.329e+05 ... 4.762e+05 4.762e+05
    spatial_ref  int32 4B 32628
Data variables:
    count_wet    (time, y, x) int16 4MB 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0
    count_clear  (time, y, x) int16 4MB 30 30 30 30 30 30 ... 26 26 26 26 26 26
    frequency    (time, y, x) float32 9MB 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
Attributes:
    crs:           epsg:32628
    grid_mapping:  spatial_ref

Tracer le résumé du WOfS

Traçons les données WOfS pour avoir une idée des objets que nous allons transformer. Dans le code ci-dessous, nous sélectionnons d’abord les pixels où le satellite a observé de l’eau au moins 25 % de l’année, afin de pouvoir isoler les plans d’eau les plus persistants et réduire une partie du bruit avant de vectoriser le raster.

[4]:

# Select pixels that are classified as water > 25 % of the year
water_bodies = ds.frequency > 0.25

# Plot the data
water_bodies.plot(size=5)

[4]:

<matplotlib.collections.QuadMesh at 0x7fca2d2a7710>
../../../_images/sandbox_notebooks_Frequently_used_code_Rasterise_vectorise_12_1.png

Vectorisation d’un « xarray.DataArray »

Pour convertir notre objet « xarray.DataArray » en un « geopandas geodataframe » basé sur un vecteur, nous pouvons utiliser la fonction DE Africa « xr_vectorize » <https://docs.digitalearthafrica.org/en/latest/sandbox/notebooks/Tools/gen/deafrica_tools.spatial.html#deafrica_tools.spatial.xr_vectorize>`__ dans le module deafrica_tools.spatial. Cet outil est basé sur la fonction rasterio.features.shape et peut accepter n’importe lequel des arguments de rasterio.features.shape en utilisant la même syntaxe.

Dans la cellule ci-dessous, nous utilisons l’argument « mask=water_bodies.values==1 » pour indiquer que nous souhaitons uniquement convertir les valeurs de l’objet xarray qui sont égales à 1.

Remarque : « xr_rasterize » et « xr_vectorize » tenteront d’obtenir automatiquement les « crs » et « transform » à partir des données d’entrée, mais si les données ne contiennent pas ces informations, vous devrez les fournir manuellement. Dans la cellule ci-dessous, nous obtiendrons les « crs » et « transform » à partir de l’ensemble de données d’origine.

[5]:

gdf = xr_vectorize(water_bodies,
                   crs=ds.crs,
                   mask=water_bodies.values==1)

print(gdf.head())

   attribute                                           geometry
0        1.0  POLYGON ((462495 1514655, 462495 1514625, 4625...
1        1.0  POLYGON ((465945 1514655, 465945 1514625, 4659...
2        1.0  POLYGON ((466125 1514655, 466125 1514625, 4661...
3        1.0  POLYGON ((466275 1514655, 466275 1514625, 4663...
4        1.0  POLYGON ((462315 1514625, 462315 1514565, 4623...

Tracer notre raster vectorisé

[6]:

gdf.plot(figsize=(6, 6))

[6]:

<Axes: >
../../../_images/sandbox_notebooks_Frequently_used_code_Rasterise_vectorise_16_1.png

Exporter en tant que fichier de formes

Notre fonction nous permet également d’exporter très facilement le « GeoDataFrame » en tant que « shapefile » pour l’utiliser dans d’autres applications en utilisant le paramètre « export_shp ».

[7]:

gdf = xr_vectorize(da=water_bodies,
                   crs=ds.crs,
                   mask=water_bodies.values == 1.,
                   output_path='test.shp')

Exporting vector data to test.shp

Rastérisation d’un shapefile

En utilisant la fonction `xr_rasterize <https://docs.digitalearthafrica.org/en/latest/sandbox/notebooks/Tools/gen/deafrica_tools.spatial.html#deafrica_tools.spatial.xr_rasterize>`__ dans le module deafrica_tools.spatial (basée sur la fonction rasterio : rasterio.features.rasterize, et peut accepter n’importe lequel des arguments dans rasterio.features.rasterize en utilisant la même syntaxe), nous pouvons transformer le geopandas.GeoDataFrame en un « xarray.DataArray ».

Comme nous avons déjà chargé le « GeoDataFrame », nous n’avons pas besoin de lire le shapefile, mais si nous voulions d’abord lire un shapefile, nous pouvons utiliser gpd.read_file().

Cette fonction utilise un objet « xarray.dataArray » comme modèle pour convertir le « geodataframe » en un objet raster (le modèle fournit la « taille », « crs », « dimensions », « transform » et les « attributs » du tableau de sortie).

[8]:

water_bodies_again = xr_rasterize(gdf=gdf,
                                  da=water_bodies,
                                  crs=ds.crs)

print(water_bodies_again)

<xarray.DataArray (y: 1478, x: 1447)> Size: 17MB
array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ...,
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0]])
Coordinates:
  * y            (y) float64 12kB 1.515e+06 1.515e+06 ... 1.47e+06 1.47e+06
  * x            (x) float64 12kB 4.328e+05 4.329e+05 ... 4.762e+05 4.762e+05
    spatial_ref  int32 4B 32628

Exporter au format GeoTIFF

xr_rasterize permet également d’exporter les résultats au format GeoTIFF en utilisant le paramètre export_tiff. Pour ce faire, un tableau named est requis. Si aucun n’est fourni, la fonction en fournira un par défaut.

[9]:

water_bodies_again = xr_rasterize(gdf=gdf,
                                  da=water_bodies,
                                  crs=ds.crs,
                                  output_path='test.tif')

Exporting raster data to test.tif

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 :

[10]:

print(datacube.__version__)

1.8.20

Dernier test :

[11]:

from datetime import datetime
datetime.today().strftime('%Y-%m-%d')

[11]:

'2025-01-15'