Suite

Créer un raster avec des valeurs de cellules de grille basées sur la zone couverte par des polygones

Créer un raster avec des valeurs de cellules de grille basées sur la zone couverte par des polygones


J'ai des données de polygone avec une superficie couverte par les forêts (les données sont ici - https://www.dropbox.com/s/zgckliydalljw6a/sp_data.zip?dl=0). Je veux convertir des polygones en raster. La valeur de chaque cellule de la grille doit être basée sur la zone couverte par le polygone. Par exemple, si la taille de la cellule de la grille est de 100 m x 100 m (10 000 m2) et que les polygones couvraient 5 000 m2 de cette cellule, la valeur devrait être 50 (ou 0,5). Ci-dessous, vous pouvez trouver ma proposition de solution:

Données brutes:

bibliothèque('rgdal') bibliothèque('raster') bibliothèque('rgeos') bibliothèque('maptools') sp_data <- readOGR(".", "sp_data") spplot(sp_data)

J'ai créé un nouveau raster avec une résolution de 100 mètres et une étendue égale à la zone de délimitation du polygone :

r <- raster() bb <- extent(290000, 300000, 500000, 510000) extent(r) <- bb res(r) <- 100

J'ai également créé un polygone d'arrière-plan - avec les valeurs des zones sans forêt égales à 0. Ensuite, j'ai fusionné ces deux données de polygone :

sp_back <- as(extent(r), "SpatialPolygons") proj4string(sp_back) <- proj4string(sp_data) sp_back <- gDifference(sp_back, sp_data) plot(sp_back, col="red") plot(sp_data, col=" green") sp_back <- SpatialPolygonsDataFrame(sp_back, data=data.frame(value=rep(0, length(sp_back)), row.names=row.names(sp_back))) sp_back <- spChFIDs(sp_back, "new_id" ) sp_bind <- spRbind(sp_data, sp_back) spplot(sp_bind)

De plus, j'ai créé un nouveau raster avec une résolution réduite (500 mètres):

r2 <- r res(r2) <- 500

J'ai pixellisé les données de polygone dans le premier raster :

sp_raster <- rasterize(sp_bind, r, field="value", fun=max) spplot(sp_raster, aspect="iso")

À la fin, j'ai rééchantillonné les premières valeurs raster dans le deuxième raster :

sp_raster2 <- resample(sp_raster, r2, method="bilinear") spplot(sp_raster2, aspect="iso")

Sur la base des résultats, j'ai quelques questions:

  1. Tout d'abord, ma solution et le résultat sont-ils corrects ?
  2. La conversion entre vecteur et raster simplifie en quelque sorte mes données (par exemple, perte de petits polygones) et ajoute probablement quelques erreurs. Le résultat du rééchantillonnage est affecté par cette simplification. Existe-t-il une solution alternative pour calculer la part de la surface du polygone dans la cellule de la grille ?

Les réponses utilisant R, Grass ou SAGA seront les bienvenues.


Puisque vous avez mis GRASS dans les balises, voici une solution basée sur GRASS :

Tout d'abord, vous devez savoir exactement dans quel système de coordonnées se trouvent les données d'origine (comme toujours avec GRASS). Je vois que le fichier *.prj contient "TRANSVERSE MERCATOR" mais ce n'est pas l'une des zones UTM standard. Puisque vous n'avez pas mentionné le code EPSG de ces données, voici la chaîne proj4 pour créer un nouvel emplacement GRASS correspondant :

+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=WGS84 +units=m +no_defs

Maintenant, démarrez GRASS, créez un LOCATION (basé sur ce qui précède) et un MAPSET.

Et maintenant, importez le fichier de formes dans GRASS en exécutant :

v.in.ogr input="sp_data.shp" output=sp_data

Ajoutez une colonne pour la surface de chaque entité et calculez sa surface :

v.db.addcolumn map=sp_data column="area_sqm INT" v.to.db sp_data option=area column=area_sqm unit=meters

Convertissez maintenant en raster, en utilisant la colonne area_sqm comme valeur raster. Mais d'abord définir la région de calcul GRASS:

g.region -p vector=sp_data res=100

Créez maintenant une grille vectorielle de polygones à la résolution raster :

v.mkgrid map=grid

Intersection avec les polygones forestiers et calcul des zones de ces polygones d'intersection : Intersectez les polygones de la forêt, dissolvez les polygones à l'intérieur de la même cellule de grille, calculez les zones de ces polygones d'intersection et déplacez la valeur de zone vers les cellules de grille :

v.overlay ain=grid bin=sp_data operator=and output=forest_grid_intersect v.dissolve --overwrite [email protected] column=a_cat output=forest_grid_intersect_dissolve v.db.addtable map=forest_grid_intersect_dissolve colonnes="va.to.sqm INT" db forest_grid_intersect_dissolve option=area column=area_sqm unit=meters v.db.join [email protected] column=cat other_table=forest_grid_intersect_dissolve other_column=cat

v.db.addcolumn map=forest_grid_intersect column="area_sqm INT" v.to.db forest_grid_intersect option=area column=area_sqm unit=meters

et convertissez en raster, en utilisant l'attribut area comme valeur raster :

v.to.rast input=grid type=area output=sp_data_rast use=attr attribute_column=area_sqm

Et tu as fini. J'espère avoir raison cette fois :-) .


J'essayais juste de comprendre comment faire la même chose en r. Voici l'approche que j'ai trouvé :

-Créez un raster qui a les dimensions et la taille de cellule que vous voulez

-Utilisez la fonction 'rasterize' du package raster, en utilisant l'option "getCover"

library(raster) zones humides <- readShapePoly("data/wetlands.shp", proj4string = CRS("+init=EPSG:26907)) r <- raster(ncols = 200, nrows = 200, resolution = 5, crs = CRS ("+init=EPSG:26907)) zones humidesAREA <- rasterize(zones humides, r, getCover = TRUE)

Je pense que cela fait essentiellement le même processus que dans le message d'origine. La fonction rasterize divise la cellule raster en 100 sous-unités, et la sortie est la proportion de celles qui étaient couvertes par le polygone.


Voir la vidéo: monikulmio