3

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.

7

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)
  • thanks a lot for your help! The Transformation works partly: extent of the shapefile is ok. CRS is obviously same as raster. BUT my extracted values are unfortunately not correct. (i do a manual check by clicking pixel values in QGIS- I know that as well could be a fault - but some of my fields have very obvious borders so I can do a nice check.) so in the end I assign another projection for extracting my values - which is the on the fly projection from QGIS (next comment) – kaltfront Oct 8 '14 at 15:20
  • 1
    (+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") with this shp file I get the correct values but my extent is that one I mentioned in the beginning. processing the command: e1<-extract(rst, shputm, fun=mean, factors=TRUE) I also get this Warning message: In .local(x, y, ...) : Transforming SpatialPolygons to the CRS of the Raster So obviously there has to be some transformation done in order to do the extraction. – kaltfront Oct 8 '14 at 15:23
  • (I dont get warning message BUT false values when the shapefile with the crs(raster) is used.) – kaltfront Oct 8 '14 at 15:24
  • Saved me hours. Very very useful. – Ritesh Porwal Aug 22 at 12:18
0

If you want to extract value from a SRTM data then you should use a point shapefile because raster image change their pixel value at different level and in polygon layer you can not find right value. so you have to need a point shapfile to extract value.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.