Suite

Reprojection du shapefile dans R et extraction finale des valeurs

Reprojection du shapefile dans R et extraction finale des valeurs


J'ai besoin d'extraire des valeurs de certains fichiers raster (pile) à l'aide d'un fichier de formes. pour les mettre tous les deux dans le même CRS, j'ai d'abord transformé le shapefile en utilisant les codes EPSG. En l'ouvrant dans un logiciel SIG (QGIS, Saga), il n'était pas au bon endroit. (également en récupérant les rasters CRS et en reprojetant le shapefile).

Ensuite, j'ai extrait manuellement le code crs correct (+proj=lcc +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.3333333333333333 +x_0=400000 +y_0=400000 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units =m +no_defs) de la "projection à la volée" dans QGis pour mon R-Code.

Au final, j'ai trouvé ce qui suit :

shp<-readOGR(dsn="hagel", layer="testflaechen_20131203") crs(shp) # obtenir le système de coordonnées de référence (CRS) >arguments CRS : +proj=lcc +lat_1=46 +lat_2=49 +lat_0=47.5 + lon_0=13.333333333333333 +x_0=400000 +y_0=400000 +ellps=bessel +units=m +no_defs shputm<-shp crs(shputm)<-"+proj=lcc +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0 =13.3333333333333 +x_0=400000 +y_0=400000 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs" compareCRS(shp, shputm) [1] TRUE writeOGR(shputm , dsn="hagel", layer="testflaechen_20131203_utm", driver="ESRI Shapefile", overwrite_layer=TRUE)

Alors la chose étrange est quecomparerCRSme dit qu'ils sont identiques. Dans QGis, ils ne semblent pas l'être et aussi si j'utiliseextrait()dans d'autres étapes avec shp et shputm cela me donne des valeurs différentes.

Un autre problème ou indication que quelque chose ne va pas avec la projection peut être vu dansextshp1<-extent(shputm)conduisant aux valeurs suivantes :

classe extshp1 : Etendue xmin : 649456,4 xmax : 657021,6 ymin : 488263,9 ymax : 505175,5

Bien que les valeurs devraient ressembler à celles-ci :extshp<-extent(c(623874.1, 634993, 5343486, 5362521))(obtenu les valeurs manuellement viadrawextent())

shapefiles sur dropbox et un exemple de raster de la région

En fait, tout cela devrait être simple, mais j'ai vraiment du mal avec ça.


Le code ci-dessus ressemble plus à un cloneshpdansshputmpuis affecter la sortie decrs(shp)àshputmsans effectuer une reprojection réelle. Quoi qu'il en soit, si vous importez à la fois le fichier de formes et le raster NDVI, puis reprojetezshputilisantspTransform, l'extraction de données ultérieure devrait fonctionner correctement. Aussi, la sortie deétendue(shp_utm)correspond à peu près à l'étendue du fichier de formes extrait manuellement que vous avez décrit ci-dessus.

# Packages obligatoires library(rgdal) library(raster) # Shapefile import shp <- readOGR(dsn = "data", layer = "testflaechen_20131203") # crs(shp) # Raster import rst <- raster("data/ndvi_LE71890272009095ASN00.tif ") # crs(rst) # reprojection du fichier de formes shp_utm <- spTransform(shp, crs(rst)) # compareCRS(rst, shp_utm) # extent(shp_utm) # Extraction de la valeur extract(rst, shp_utm)

Si vous souhaitez extraire la valeur d'une donnée SRTM, vous devez utiliser un fichier de formes de points car l'image raster modifie sa valeur de pixel à différents niveaux et dans la couche de polygones, vous ne pouvez pas trouver la bonne valeur. vous devez donc avoir besoin d'un fichier shapfile pour extraire la valeur.


Voir la vidéo: R studio: extract raster values by XY coordinate and line