【问题标题】:Raster not cropping to shapefile栅格未裁剪为 shapefile
【发布时间】: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


【解决方案1】:
library(rgdal)
library(raster)
library(rgeos)

usa.states <- readOGR("states.shp", layer = "states")

# Here we subset once
Corners <- usa.states[usa.states$STATE_NAME %in% c("Colorado", "New Mexico","Arizona","Utah"),]

# Dissolve polygons into one
Corners <- gUnaryUnion(Corners)

# Create a 20x20 raster using the extent of Corners
# The number of rows and columns can be change to increase/reduce the resolution
r <- raster(extent(Corners), ncol=20, nrow=20, crs=CRS(proj4string(Corners)))

# Rasterize
Corners.crop <- rasterize(Corners, r)

【讨论】:

  • 大部分情况下的答案都很好,但是我真的需要 1Km 的最终分辨率。有什么方法可以准确地设置它(使纬度线的不等间距生效)还是更容易近似和稍后更改?
  • 如所写,“可以更改行数和列数以增加/减少分辨率”。有了这个预测,它将在 1 公里左右。我让你根据经纬度的范围来计算行数和列数。
猜你喜欢
  • 2018-09-15
  • 1970-01-01
  • 2014-05-29
  • 2022-09-24
  • 2019-01-06
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多