I have a series of binary raster layers (ascii file) showing presence/absence of a species in Europe and Africa. The file is based on unprojected lat/long (WGS84) data. My aim is to calculate the area of presence using R (I don't have access to ArcGIS).
I know that the raster package has a function for calculating area, but I'm worried that this won't be accurate for unprojected data. I have also looked at the cellStats function in the raster package, and can use this to "sum" the number of cells occupied, but I feel this has the same problem.
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
Anybody know of a way to calculate the presence area accurately, accounting for the earth curvature?
Thanks!
I do not know what is going in with your file (why you get warning #2). But is here is a work around
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_Robert\2015-01-26_213612_1208_85354.grd
names : layer
values : 0.2227891, 0.8605576 (min, max)
Now you have the area of each cell in km2. By multiplying these values with Raster objects with presence/absence values and then using cellStats( , 'sum') you can obtain the total area with presence.