使用该raster
包,我读入了两个数据集,一个 ASCII 栅格和一个 ESRI shapefile。我想将栅格(水温)数据提取到 shapefile 的全部范围,即湖岸线。
使用该函数SpatialPolygonsDataFrame
读入时,ESRI shapefile 被视为一个。shapefile()
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
栅格的(即使用crs
函数中的参数强制转换为以下内容raster()
):
+proj=merc +datum=WGS84 +ellps=WGS84 +units=m +towgs84=0,0,0
然后我使用包中的spTransform()
函数rgdal
将 shapefile 的空间参考强制转换为栅格的空间参考。
spTransform(shapefile, CRS(projection(raster)))
最后,我提交了以下内容:
extract(raster,shapefile,method="simple",fun=mean,small=TRUE,na.rm=TRUE,df=FALSE)
但是,extract()
只是返回类型NULL
的对象list
。我认为这个问题是由坐标引用的显式强制引起的。
此外,以下是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
我在这些论坛上搜索了大量类似的问题,但没有得到解决。有人可以借我一个(虚拟的)手吗?
非常感谢,提前。