Map raw data and mean data based on the shapefile

Go To


sI 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
names(pts) <- c( "x","y","data")

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

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")


Now I would like to create a separate map which shows the means of pts across different regions. The shapefile I want to use is from , however, this is a north america map. How can I create the map showing only US based on this north america map? Or is there another better way to do this? thanks so much.

2012-04-04 19:39
by Dan


I think that cutting out the non-US data based on the data in the shapefile alone would be hard, since the regions do not correspond to political boundaries - that could be done with rgeos though.

Assuming that "eco" is a SpatialPolygonsDataFrame read in by rgdal::readOGR or maptools::readShapeSpatial, see the available key data for indexing:

sapply(, function(x) if(!is.numeric(x)) unique(x) else NULL)

If you just want to plot it, set up a map with only the US region to start with and then overplot.

map("usa", col = "transparent")

We see that the data is in Lambert Azimuthal Equal Area:

[1] " +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"

So require(rgdal) eco.laea <- spTransform(eco, CRS("+proj=longlat +ellpse=WGS84")) plot(eco.laea, add = TRUE)

If you want to plot in the original Lambert Azimuthal Equal Area you'll need to get the bounding box in that projection and start the plot based on that, I just used existing data to make an easy example. I'm pretty sure the data could also be cropped with rgeos against another boundary too, but depends what you actually want.

2012-04-04 22:12
by mdsumner
Thanks for your input. I get what you mean by doing this. However, how would you plot the means on the region map - Dan 2012-04-04 23:24
Use over() in sp to identify regions for each point, then aggregate the values. Your question does not actually ask for that explicitly - mdsumner 2012-04-04 23:35
Sorry about not making my question explicitly. Could you please be more specific? thank - Dan 2012-04-05 14:38