从二进制未投影栅格计算占用面积
Calculating area of occupancy from a binary unprojected raster
我有一系列二进制栅格图层(ascii 文件)显示 presence/absence 欧洲和非洲的一个物种。该文件基于未投影的 lat/long (WGS84) 数据。我的目标是使用 R 计算存在区域(我无权访问 ArcGIS)。
我知道栅格包有计算面积的功能,但我担心这对未投影的数据不准确。我也看过raster包中的cellStats函数,可以用它来"sum"占用的像元数,但我觉得这也有同样的问题。
jan<-raster("/filelocation/file.asc")
jan
class : RasterLayer
dimensions : 13800, 9600, 132480000 (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -20, 60, -40, 75 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : "/filelocation"
names : file.asc
values : -2147483648, 2147483647 (min, max)
area(jan)
class : RasterLayer
dimensions : 13800, 9600, 132480000 (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -20, 60, -40, 75 (xmin, xmax, ymin, ymax)
coord. ref. : NA
names : layer
values : 6.944444e-05, 6.944444e-05 (min, max)
Warning messages:
1: In .local(x, ...) :
This function is only useful for Raster* objects with a longitude/latitude coordinates
2: In .rasterFromRasterFile(grdfile, band = band, objecttype, ...) :
size of values file does not match the number of cells (given the data type)
cellStats(jan,"sum")
[1] 3559779
任何人都知道一种准确计算存在区域的方法,考虑到地球曲率?
谢谢!
我不知道您的文件中发生了什么(为什么您会收到警告 #2)。但是这里有一个变通方法
r <- raster(nrow=13800, ncol=9600, xmn=-20, xmx=60, ymn=-40, ymx=75)
# equivalent to r <- raster(jan)
x = area(r)
x
class : RasterLayer
dimensions : 13800, 9600, 132480000 (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -20, 60, -40, 75 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84
data source : c:\temp\R_raster_Robert15-01-26_213612_1208_85354.grd
names : layer
values : 0.2227891, 0.8605576 (min, max)
现在您得到了以 km2 为单位的每个单元格的面积。通过将这些值与具有 presence/absence 值的 Raster 对象相乘,然后使用 cellStats( , 'sum') 您可以获得存在的总面积。
我有一系列二进制栅格图层(ascii 文件)显示 presence/absence 欧洲和非洲的一个物种。该文件基于未投影的 lat/long (WGS84) 数据。我的目标是使用 R 计算存在区域(我无权访问 ArcGIS)。
我知道栅格包有计算面积的功能,但我担心这对未投影的数据不准确。我也看过raster包中的cellStats函数,可以用它来"sum"占用的像元数,但我觉得这也有同样的问题。
jan<-raster("/filelocation/file.asc")
jan
class : RasterLayer
dimensions : 13800, 9600, 132480000 (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -20, 60, -40, 75 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : "/filelocation"
names : file.asc
values : -2147483648, 2147483647 (min, max)
area(jan)
class : RasterLayer
dimensions : 13800, 9600, 132480000 (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -20, 60, -40, 75 (xmin, xmax, ymin, ymax)
coord. ref. : NA
names : layer
values : 6.944444e-05, 6.944444e-05 (min, max)
Warning messages:
1: In .local(x, ...) :
This function is only useful for Raster* objects with a longitude/latitude coordinates
2: In .rasterFromRasterFile(grdfile, band = band, objecttype, ...) :
size of values file does not match the number of cells (given the data type)
cellStats(jan,"sum")
[1] 3559779
任何人都知道一种准确计算存在区域的方法,考虑到地球曲率?
谢谢!
我不知道您的文件中发生了什么(为什么您会收到警告 #2)。但是这里有一个变通方法
r <- raster(nrow=13800, ncol=9600, xmn=-20, xmx=60, ymn=-40, ymx=75)
# equivalent to r <- raster(jan)
x = area(r)
x
class : RasterLayer
dimensions : 13800, 9600, 132480000 (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -20, 60, -40, 75 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84
data source : c:\temp\R_raster_Robert15-01-26_213612_1208_85354.grd
names : layer
values : 0.2227891, 0.8605576 (min, max)
现在您得到了以 km2 为单位的每个单元格的面积。通过将这些值与具有 presence/absence 值的 Raster 对象相乘,然后使用 cellStats( , 'sum') 您可以获得存在的总面积。