Coordinate Reference System of Raster & Polygons in R

Question

I have several polygons & I like to extract mean values from several raster layers within these polygons. When I added those to ArcMap I realized that the projections of the two data types do not match. I could solve the problem for the display in ArcGIS by using the Project tool (Data Management toolbox > Projections and Transformations toolset > Raster). So I tried to standardize the projection by loading the data into R in the following manner (part of the code):

Rasters:

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 plot the polygons and the rasters seperately, but when I try to plot one of the polygons over a 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 still do not "overlap". This is also indicated by the different bounding boxes, I think:

> 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

Because I want to extract the raster values within the polygons, I have to be sure that they are referenced in a correct way. Does someone have an idea what the problem could be?

Answer 1

Your polygon shapefile has a coordinate system that isnt lat-long - the numbers are very large and probably metres in some system. Assigning a proj4string won't reproject the data to lat-long, it just sets the label of what coordinates it thinks it is. In this case, its wrong!

You need to make sure your polygons get the right proj4string for the numbers they have in them - there may be a [shapefile].prj file along with the .shp and .dbf that tells you. Set the proj4string to that.

Then you can use spTransform from sp or rgdal to project your polygons to lat-long WGS84 coordinates.

Its always best to transform polygons to raster coordinates, since messing with raster coordinates can mean reprojecting the whole grid, which is messy...