I have several polygons, and I like to extract averages from several raster layers inside these polygons. When I added them to ArcMap, I realized that the predictions of the two data types do not match. I could solve the problem for display in ArcGIS using the Project tool (Data Management Tool> Projections and Transform Tools> Raster). So I tried to standardize the projection by loading the data into R as follows (part of the code):
raster images on:
for (i in 1:length(rasterlist1)) {ndvi_raster_stack1[i]<-raster(rasterlist1[i]) raster::NAvalue(ndvi_raster_stack1[[i]])<--999 projection(ndvi_raster_stack1[[i]])<-"+proj=utm +ellps=WGS84 +datum=WGS84 +units=m"} > ndvi_raster_stack1[[1]] class : RasterLayer dimensions : 226, 150, 33900 (nrow, ncol, ncell) resolution : 0.57504, 0.5753628 (x, y) extent : -28.728, 57.528, -55.08, 74.952 (xmin, xmax, ymin, ymax) coord. ref. : +proj=utm +ellps=WGS84 +datum=WGS84 +units=m +towgs84=0,0,0 values : Z:\master\lusmeg_sw_kernel_data\ndvi0910\Y2008_P47.tif min value : -91 max value : 550.8125
Polygons:
for (i in 1:length(poplist)) {pop_kernels[i]<-readShapeSpatial(poplist[i],repair=TRUE,proj4string=CRS("+proj=utm +ellps=WGS84 +datum=WGS84 +units=m")) pop_kernels[[i]]<-unionSpatialPolygons(pop_kernels[[i]],ID=c(rep(1,times=length(pop_kernels[[i]])-1),0),threshold=NULL,avoidGEOS=FALSE)} > str(pop_kernels[[1]]) Formal class 'SpatialPolygons' [package "sp"] with 4 slots ..@ polygons :List of 2 .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots .. .. .. ..@ Polygons :List of 2 .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots .. .. .. .. .. .. ..@ labpt : num [1:2] 2404422 893343 .. .. .. .. .. .. ..@ area : num 1.15e+12 .. .. .. .. .. .. ..@ hole : logi FALSE .. .. .. .. .. .. ..@ ringDir: int 1 .. .. .. .. .. .. ..@ coords : num [1:1625, 1:2] 2551236 2533236 2533236 2523236 2523236 ... .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2 .. .. .. .. .. .. .. .. ..$ : NULL .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y" .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots .. .. .. .. .. .. ..@ labpt : num [1:2] 2468549 865776 .. .. .. .. .. .. ..@ area : num 6.31e+11 .. .. .. .. .. .. ..@ hole : logi TRUE .. .. .. .. .. .. ..@ ringDir: int -1 .. .. .. .. .. .. ..@ coords : num [1:1385, 1:2] 2551236 2551236 2563236 2563236 2569236 ... .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2 .. .. .. .. .. .. .. .. ..$ : NULL .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y" .. .. .. ..@ plotOrder: int [1:2] 1 2 .. .. .. ..@ labpt : num [1:2] 2404422 893343 .. .. .. ..@ ID : chr "0" .. .. .. ..@ area : num 1.15e+12 .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots .. .. .. ..@ Polygons :List of 1 .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots .. .. .. .. .. .. ..@ labpt : num [1:2] 2468549 865776 .. .. .. .. .. .. ..@ area : num 6.31e+11 .. .. .. .. .. .. ..@ hole : logi FALSE .. .. .. .. .. .. ..@ ringDir: int 1 .. .. .. .. .. .. ..@ coords : num [1:1385, 1:2] 2551236 2541236 2541236 2529236 2529236 ... .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2 .. .. .. .. .. .. .. .. ..$ : NULL .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y" .. .. .. ..@ plotOrder: int 1 .. .. .. ..@ labpt : num [1:2] 2468549 865776 .. .. .. ..@ ID : chr "1" .. .. .. ..@ area : num 6.31e+11 ..@ plotOrder : int [1:2] 1 2 ..@ bbox : num [1:2, 1:2] 1819236 207017 3013236 1577017 .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : chr [1:2] "x" "y" .. .. ..$ : chr [1:2] "min" "max" ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots .. .. ..@ projargs: chr " +proj=utm +ellps=WGS84 +datum=WGS84 +units=m +towgs84=0,0,0"
I can split polygons and rasters separately, but when I try to build one of the polygons above the raster, the polygons are not displayed:
plot(ndvi_raster_stack1[[1]],xlab="Longitude",ylab="Latitude") plot(pop_kernels[[1]],col="black",add=TRUE)
It seems that they are still not "overlapping." This is also evidenced by various limiting fields:
> bbox(ndvi_raster_stack1[[1]]) min max s1 -28.728 57.528 s2 -55.080 74.952 > bbox(pop_kernels[[1]]) min max x 1819236 3013236 y 207017 1577017
Since I want to extract raster values in polygons, I need to be sure that they are referenced correctly. Can anyone understand what the problem is?