1

I have the dataset (pts) like this:

x <- seq(-124.25,length=115,by=0.5)    
y <- seq(26.25,length=46,by=0.5)
z = 1:5290

longlat <- expand.grid(x = x, y = y)  # Create an X,Y grid
pts=data.frame(longlat,z) 
names(pts) <- c( "x","y","data")

I knew that I can map the dataframe (pts） into a map by doing:

library(sp)
library(rgdal)
library(raster)
library(maps)
coordinates(pts)=~x+y

proj4string(pts)=CRS("+init=epsg:4326") # set it to long, lat

pts = spTransform(pts,CRS(" +init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
pts <- as(pts, "SpatialPixelsDataFrame")
r = raster(pts)
projection(r) = CRS(" +init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")

plot(r)
map("usa",add=T)

my question is that how I can calculate the means of the 10 EPA Regions and map the means?

the EPA regions can be found at the bottom of the webpage at http://www.epa.gov/wed/pages/ecoregions/level_iii_iv.htm

Thanks

Dan
  • 435
  • 1
  • 8
  • 14
  • Before even getting to your question, there are several problems with your code. I've fixed one by adding a call to `library(raster)`. Then, if you want `r` to contain the values in `z`, you'll need to also add the following line before your call to `raster(pts)`: pts <- as(pts, "SpatialPixelsDataFrame"). (Because you haven't done this, `plot(r)` currently fails with an error message). Finally, to make this an easily reproducible example, you could also add code to create a simple SpatialPolygon object that overlays the RasterLayer. (Otherwise, folks have to create one in order to help you). – Josh O'Brien Mar 30 '12 at 00:17
  • I can't get rgdal to work with R 2.14.2. it says there's a name space issue, that indicates either I'm messing up or the package hasn't been updated and you're using an older version of R. So the best help I can give is a link to a question I asked that may be similar [(LINK)](http://stackoverflow.com/questions/9441778/improve-centering-county-names-ggplot-maps) – Tyler Rinker Mar 30 '12 at 00:17
  • @TylerRinker -- I'm also running R 2.14.2, on Windows XP, as well as the latest version of **rgdal** ('0.7.8', published 2012-01-18), and it works without a hitch. Out of curiosity, are you on Windows, or something else? – Josh O'Brien Mar 30 '12 at 00:25
  • Win 7 but I've been having problems with my computer lately. That may be it. Little glitches. I'm shopping for a new one (8GB memory this time) – Tyler Rinker Mar 30 '12 at 00:48
  • Josh, I add pts <- as(pts, "SpatialPixelsDataFrame") in my code. Thanks for catching that. – Dan Mar 30 '12 at 01:10
  • I know you want the means, but are you sure? The simple means will not calculate the center of the polygon. It will calculate the mean position of each edge. If you have many edges in one side (for example, in the coastline) and not as many in another side (as in Texas), this would skew your center significantly. – dmvianna Dec 06 '12 at 23:39

1 Answers1

2

First read the shape file

er <- readOGR("Eco_Level_III_US.shp", "Eco_Level_III_US")

Then make sure ther raster r and the ecoregions er have the same projection

er.4326 <- spTransform(er, CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))

Extract the r raster data from the shapefile (that may take a few minutes), calculate the mean and add it to your polyongs.

er.v <- extract(r, er.4326)
means <- sapply(er.v, mean)
er.4326$means <- means

And finally plot it

 spplot(er.4326, "means")
johannes
  • 14,043
  • 5
  • 40
  • 51
  • Awesome. how would it be possible to provide a legend for each region which is correspondent to the ones shown on their original map? Sorry about the lame questions. – Dan Mar 30 '12 at 14:51
  • I noticed that there are 31 regions on this map. Is there a way that we can calculate the means of the polygons shown on the map at the bottom their webpage? I could not find a shapefile for that map. – Dan Mar 30 '12 at 15:00