【发布时间】:2015-02-26 05:54:28
【问题描述】:
我正在尝试合并来自shapefile 的一些状态,并生成一个我可以使用downstream 的raster。我已经合并了状态,但是当我创建一个空栅格以使用裁剪功能进行栅格化时似乎失败了。我对R 中的GIS features 很陌生,非常感谢您的帮助。
Shapefile is fromhttp://www.arcgis.com/home/item.html?id=f7f805eb65eb4ab787a0a3e1116ca7e5
library(maptools)
library(shapefiles)
library(raster)
usa.states <- readOGR(dsn = "states_21basic/", layer = "states")
head(usa.states)
Co=usa.states[usa.states@data$STATE_NAME== "Colorado",]
Nm=usa.states[usa.states@data$STATE_NAME== "New Mexico",]
Az=usa.states[usa.states@data$STATE_NAME== "Arizona",]
Ut=usa.states[usa.states@data$STATE_NAME== "Utah",]
Corners= spRbind(spRbind(spRbind(Co,Ut),Nm),Az)
CRS="+proj=longlat +datum=WGS84"
Corners=spTransform(Corners, CRS(CRS))
> extent(Corners)
class : Extent
xmin : -114.8218
xmax : -102.0372
ymin : 31.33563
ymax : 42.0023
cor.ext=extent(Corners)
r<-raster(ncol=ncol(Corners), nrow=nrow(Corners), crs=CRS)
Corners.crop= crop(r,cor.ext, snap="out")
当我随后调用'Corners.crop' 的范围时,我收到了:
> extent(Corners.crop)
class : Extent
xmin : -180
xmax : -36
ymin : 0
ymax : 45
我很困惑我缺少什么才能让它工作。 我也希望有一个 1Km 的分辨率,并且很好奇在空栅格上或在我栅格化形状之后更改分辨率是否会更好。
【问题讨论】:
-
readOGR在rgdal包中,crop在raster包中。无需重新投影数据,它们已经在+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0中。你检查ncol(Corners)和nrow(Corners)的输出了吗? -
@Pascal 确实遗漏了我安装的一些软件包,对不起。我知道我真的不需要重新投影数据,但对于最终产品来说,这似乎是一个容易的地方。就“ncol”和“nrow”而言,该部分将被替换为可以让我获得更好分辨率的内容。
标签: r raster shapefile maptools