[GIS] Reprojection of shapefile in R and final extraction of values

coordinate systemrrastershapefile

I need to extract values from some raster files (stack) using a shapefile. to bring both of them in the same CRS, I first transformed the shapefile using the EPSG codes. Opening it in an GIS software (QGIS, Saga) it was not in the correct place. (also via retrieving the rasters CRS and reprojecting the shapefile).

Afterwards I extracted manually the correct crs code (+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) from the "on the fly projection" in QGis for my R-Code.

In the end I came up with the following: 

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)

So the strange thing is that compareCRS tells me that they are the same. In QGis they don't appear to be and also if I use extract() in further steps with shp and shputm it gives me different values.

Another problem or indication that something must be wrong with the projection can be seen in extshp1<-extent(shputm) leading to following values:

extshp1
class : Extent
xmin : 649456.4
xmax : 657021.6
ymin : 488263.9
ymax : 505175.5

Although the values should be something like these: extshp<-extent(c(623874.1, 634993, 5343486, 5362521)) (got values manually via drawextent())

shapefiles on dropbox and a sample raster of the region

Actually this all should be simple, but I am really struggling with it.

Best Answer

The above code looks more like you just clone shp into shputm and then assign the output of crs(shp) to shputm without performing an actual reprojection. Anyway, if you import both the shapefile and the NDVI raster, and then reproject shp using spTransform, subsequent data extraction should work out fine. Also, the output of extent(shp_utm) roughly agrees with the manually extracted shapefile extent you depicted above. 

# 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)
Related Question