R - 'raster'包 - RasterLayer（“raster”）不提取到SpatialPolygonsDataFrame（“sp”）;只需返回NULL对象-Java 学习之路

使用 raster 包，我读入了两个数据集，一个ASCII栅格和一个ESRI shapefile . 我想将光栅的（水温）数据提取到shapefile的全部范围，shapefile是一个湖岸线 .

使用 shapefile() 函数读入时，ESRI shapefile被视为 SpatialPolygonsDataFrame .

shapefile <- shapefile("shore.shp",verbose=TRUE)

我使用 raster() 函数读入ASCII栅格 .

raster <- raster("1995_001.asc",native=TRUE,crs="+proj=merc +datum=WGS84 +ellps=WGS84 +units=m")

shapefile的坐标参考信息是：

+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0

栅格的那个（即使用 raster() 函数中的 crs 参数强制执行以下操作）：

+proj=merc +datum=WGS84 +ellps=WGS84 +units=m +towgs84=0,0,0

然后，我使用 rgdal 包中的 spTransform() 函数将shapefile的空间参考强制转换为栅格的空间参考 .

spTransform(shapefile, CRS(projection(raster)))

最后，我提交了以下内容：

extract(raster,shapefile,method="simple",fun=mean,small=TRUE,na.rm=TRUE,df=FALSE)

但是， extract() 只返回 list 类型的 NULL 对象 . 我认为这个问题是由于坐标参考的明确强制而产生的 .

另外，以下是在每个数据集上使用 show() 函数的结果：

> show(raster) class : RasterLayer dimensions : 1024, 1024, 1048576 (nrow, ncol, ncell) resolution : 1800, 1800 (x, y) extent : -10288022, -8444822, 4675974, 6519174 (xmin, xmax, ymin, ymax) coord. ref. : NA data source : in memory names : layer values : -9999, 8.97 (min, max)

> show(shapefile) class : SpatialPolygonsDataFrame features : 1 extent : 597568.5, 998261.6, 278635.3, 668182.2 (xmin, xmax, ymin, ymax) coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 variables : 3 names : AREA, PERIMETER, HECTARES min values : 59682523455.695, 5543510.075, 5968252.346 max values : 59682523455.695, 5543510.075, 5968252.346

我在这些论坛上搜索了大量类似的问题而没有解决方案 . 有人借给我（虚拟）手吗？

非常感谢，提前 .

1 回答

1
这（在空间上看似不重叠的数据）是源于坐标参考系统混淆的常见问题 . 这可能由于一些可能的原因而发生 .

1）区域实际上不重叠

2）也许你对RasterLayer的crs的赋值是错误的（或多边形的） . 请出示对象（ show(raster) ）;在询问有关 raster 包的问题时，这几乎总是有用的 .

3）您没有分配spTransform的结果 . 请记住，R通常不会进行'in-place'修改 . 一定要做 shapefile <- spTransform(shapefile, crs(raster))

始终做一个健全测试，并继续努力，直到工作正常 . 从这里开始的地方就是这样做
```
plot(raster)
plot(shapefile, add=TRUE)
```
要查看是否有任何多边形绘制在栅格数据之上 .

如果这不起作用，请查看范围 . 例如像这样（顺便说一下，这也显示了如何创建一个独立的可重现的示例/问题，其他人可以实际回答）：
```
library(raster)
# the spatial parameters of your raster
r <- raster(ncol=100, nrow=100, xmn=-10288022, xmx=-8444822, ymn=4675974, ymx=6519174, crs="+proj=merc +datum=WGS84")
values(r) <- 1:ncell(r)

# the extent of your SpatialPolygons
ep <- extent(597568.5, 998261.6, 278635.3, 668182.2)
p <- as(ep, 'SpatialPolygons')
crs(p) <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83"

# tranform both data set to longlat
rgeo <- projectRaster(r, crs='+proj=longlat +datum=WGS84')
pgeo <- spTransform(p, CRS('+proj=longlat +datum=WGS84'))

# plot on top of Google map
library(dismo)
e <- union(extent(rgeo), extent(pgeo))
g <- gmap(e,lonlat=TRUE)
plot(g, inter=TRUE)
plot(extent(rgeo), add=TRUE, col='red', lwd=2)
plot(pgeo, add=TRUE, col='blue', lwd=2)
```
显然，这两个数据源不重叠 . 只有你知道这两个中的哪一个可能是正确的，你现在可以开始研究为什么另一个在错误的地方 .
回复于 2024-05-16T08:39:49+08:00

R - 'raster'包 - RasterLayer（“raster”）不提取到SpatialPolygonsDataFrame（“sp”）;只需返回NULL对象

1 回答

相关问题