Calculating area of occupancy from a binary unprojected raster

Question

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!

It looks like you should _first_ project your data using an equal-area projection, then use `cellStats` (for example) to count the presence/absence. This question treats how to do such a projection: http://stackoverflow.com/questions/12725458/r-mapproj-lambert-azimuthal-equal-area-projection . Does it help you? — Jealie, Jan 22 '15 at 18:22
Hi @Jealie, thanks for your suggestion. However if at all possible I'd like to avoid re-projecting my data, as I have a large number of raster layers and ascii files which are all currently in the same format. — Heather, Jan 22 '15 at 19:25

score 0 · Accepted Answer · answered Jan 27 '15 at 05:40

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.

Calculating area of occupancy from a binary unprojected raster

1 Answers1