Reproyección del shapefile en R y extracción final de los valores

Question

Necesito extraer valores de unos archivos raster (stack) utilizando un shapefile. para traer ambos en el mismo CRS, primero transformé el shapefile utilizando los códigos EPSG. Al abrirlo en un software GIS (QGIS, Saga) no estaba en el lugar correcto. (también mediante la recuperación del CRS de los rasters y la reproyección del shapefile).

Después he extraído manualmente el código crs correcto ( +proj=lcc +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.33333333333333 +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 "proyección sobre la marcha" en QGis para mi R-Code.

Al final se me ocurrió lo siguiente:

shp<-readOGR(dsn="hagel", layer="testflaechen_20131203")
crs(shp) # get coordinate reference system (CRS)
>CRS arguments:
+proj=lcc +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.33333333333333 +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.33333333333333     +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)

Así que lo extraño es que compareCRS me dice que son lo mismo. En QGis no parecen serlo y además si uso extract() en otros pasos con shp y shputm me da diferentes valores.

Otro problema o indicio de que algo debe estar mal en la proyección puede verse en extshp1<-extent(shputm) lo que lleva a los siguientes valores:

extshp1 clase : Extensión xmin : 649456.4 xmax : 657021.6 ymin : 488263.9 ymax : 505175.5

Aunque los valores deberían ser algo así: extshp<-extent(c(623874.1, 634993, 5343486, 5362521)) (se obtienen los valores manualmente a través de drawextent() )

shapefiles en dropbox y un raster de muestra de la región

En realidad, todo esto debería ser sencillo, pero me está costando mucho trabajo.

Answer 1

El código anterior parece más bien que acaba de clonar shp en shputm y luego asignar la salida de crs(shp) a shputm sin realizar una reproyección real. De todos modos, si importa tanto el shapefile como el raster NDVI, y luego reproyecta shp utilizando spTransform La extracción de datos posterior debería funcionar bien. Además, la salida de extent(shp_utm) coincide más o menos con la extensión del shapefile extraído manualmente que ha representado anteriormente.

# Required packages
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)

# Shapefile reprojection
shp_utm <- spTransform(shp, crs(rst))

# compareCRS(rst, shp_utm)
# extent(shp_utm)

# Value extraction
extract(rst, shp_utm)

Answer 2

Si quiere extraer el valor de los datos SRTM, debe utilizar un archivo shapefile de puntos, ya que las imágenes rasterizadas cambian su valor de píxel en diferentes niveles y en la capa poligonal no se puede encontrar el valor correcto. Por lo tanto, se necesita un archivo shapefile de puntos para extraer el valor.