Extracting values from raster based on a shapefile

Question

I have a raster.stack of 8 environmental variables and a shapefile of Ecoregions for the southern U.S.

files<-list.files(path='E:/Ecoregions Models/Border Bioclim/',pattern='asc', full.names=TRUE)
predictors<-stack(files)

##subset the main shape file into 12 individual shapefiles for the regions of interest 
ER_11.1<-Level.2.ecoregs[Level.2.ecoregs$NA_L2CODE=="11.1",]...

I am attempting to use the extract function to pull out 5,000 random points from every ecoregion and save them as different objects.

I am having trouble projecting the individual ecoregion objects on to the raster files. I have already changed the projection but there is still a problem.

>bio_1
class       : RasterLayer 
dimensions  : 324, 444, 143856  (nrow, ncol, ncell)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -125, -88, 18.5, 45.5  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=48 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : E:\Ecoregions Models\Border Bioclim\ModelVariables\bio_1.asc 
names       : bio_1 

> ER_11.1
class       : SpatialPolygonsDataFrame 
nfeatures   : 22 
extent      : -1989641, -1450736, -1748693, -208349.6  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=48 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
nvariables  : 8
names       : NA_L2CODE,                NA_L2NAME, NA_L1CODE,                NA_L1NAME,                       NA_L2KEY,                     NA_L1KEY,   Shape_Leng,   Shape_Area 
min values  :      11.1, MEDITERRANEAN CALIFORNIA,        11, MEDITERRANEAN CALIFORNIA, 11.1  MEDITERRANEAN CALIFORNIA, 11  MEDITERRANEAN CALIFORNIA,     331.1174, 1.483274e+08 
max values  :      11.1, MEDITERRANEAN CALIFORNIA,        11, MEDITERRANEAN CALIFORNIA, 11.1  MEDITERRANEAN CALIFORNIA, 11  MEDITERRANEAN CALIFORNIA, 5417939.1114, 8.472524e+04 
>

The problem here must be in the different extents. Does anyone have any suggestions on how to resolve this in the raster package?

Robert Hijmans Robert Hijmans · Accepted Answer · 2013-07-04T20:20:09

The CRS for bio_1 is clearly wrong, You have

resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -125, -88, 18.5, 45.5  (xmin, xmax, ymin, ymax)

That is, longitude/latitude coordinates for most of the USA, but you are using

coord. ref. : +proj=utm +zone=48 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

Probably because you changed the projection name whereas you wanted to transform (project) the data? Here is what you can do:

restore the CRS of bio_1 to its original value (I presume)

projection(bio_1) <- "+proj=longlat +datum=WGS84"

transform the polygons to that same CRS

library(rgdal)
ER <- spTransform(ER_11.1, CRS(projection(bio_1)) )

And now extract the values:

v <- extract(bio_1, ER)

Extracting values from raster based on a shapefile

1 Answers