sI有这样的数据集(pts):
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")
我知道我可以通过以下操作将数据帧(pts)映射到映射中:
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)
现在我想创建一个单独的地图,显示不同地区的pts平均值。我要使用的形状文件来自ftp://ftp.epa.gov/wed/ecoregions/cec_na/NA_CEC_Eco_Level2.zip然而,这是一张北美地图。如何根据这张北美地图创建只显示美国的地图?或者还有其他更好的方法吗?非常感谢。
我认为,仅根据shapefile中的数据来剪切非美国数据是很困难的,因为这些区域不符合政治边界——不过rgeos
可以做到这一点。
假设"eco"是rgdal::readOGR
或maptools::readShapeSpatial
读入的SpatialPolygonsDataFrame
,请参阅可用的索引关键数据:
sapply(as.data.frame(eco), function(x) if(!is.numeric(x)) unique(x) else NULL)
如果你只想绘制它,可以先绘制一张只包含美国地区的地图,然后再过度绘制。
library(maps)
map("usa", col = "transparent")
我们看到数据在Lambert方位角等面积:中
proj4string(eco)
[1] " +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"
所以require(rgdal)eco.laea<-spTransform(eco,CRS("+proj=longlat+ellpse=WGS84")绘图(eco.laea,添加=真)
如果你想在原始的Lambert Azimuthal Equal Area中绘制,你需要在该投影中获得边界框,并在此基础上开始绘制,我只是使用现有数据来做一个简单的例子。我很确定数据也可以用rgeos在另一个边界上裁剪,但这取决于你实际想要什么。