Shapefile in R – How to Calculate Area of Shapefile Polygons in Square Kilometers

arearshapefile

Firstly, my shapefile has an undefined projection:

library(raster)
x <- shapefile('file.shp')

Warning message:
In .local(x, ...) : .prj file is missing

crs(x)
CRS arguments: NA

I know that my coordinates needed to be in latitude/longitude or have a distance unit in meters. For now, I still went ahead and applied my area calculation:

x$area <- area(x)/1000000

But when I calculated the total area, I got a concerning result:

sum(x$area)
[1] 3236.118

The reason for my concern is that I know for a fact that my shapefile has a total area of approximately 300 km^2

I therefore decided to define a projection using a recommended espg for where my shapefile is located (espg = 2314).

crs(x) <- CRS("+init=epsg:2314")

crs(x)
CRS arguments:
 +init=epsg:2314 +proj=cass +lat_0=10.44166666666667
+lon_0=-61.33333333333334 +x_0=86501.46392052001
+y_0=65379.0134283 +a=6378293.645208759 +b=6356617.987679838
+towgs84=-61.702,284.488,472.052,0,0,0,0 +to_meter=0.3047972654
+no_defs

However, when I reapplied my area calculation I am getting the same result as before. What am I missing?

Edit:
I'm sorry for the lack of details in the question. When I run the function print(x) I get the following properties for the shapefile (note I removed the min values and max values):

print(x)
class       : SpatialPolygonsDataFrame 
features    : 67 
extent      : 66065.81, 183871.1, 74657.79, 151865.5  (xmin, xmax, ymin, 
ymax)
crs         : NA 
variables   : 16
names       : ORDER, AREA, PERIMETER, TGOLOC_, TGOLOC_ID, NAME,     
COMM_NAME, NO_HH, MALES, FEMALES, TOTPOP, NO_BLG, NO_DWELL, NO_BUS, 
NO_INSTIT, ...

One of the attributes was indeed labelled AREA, as shown, but all values for this attribute was NA, hence the reason why I decided to try to calculate the areas myself.

Best Answer

An error in a factor of about ten in area corresponds to a scale error of about 3.16 in length (square root of ten), which I have seen happen when people use imperial foot units instead of metres (or yards for really old projection systems).

Oh look... https://epsg.io/2314

Remarks: Foot version of Trinidad 1903 / Trinidad Grid (code 30200) used by some US-based companies including Amoco Trinidad.

And your expanded EPSG string even has:

+to_meter=0.3047972654

(kicking myself for not noticing earlier...)

the projection system is in feet is not in square metres, so the area is in square feet. Scale it by feet^2 to meters^2. Your area is:

> 3236.118 * .3047972654^2
[1] 300.6398

Ta-da!

A bit more digging shows that the unit is "Clarke's Foot" which is a bit different to the International Foot, the US Surveyor's Foot, and the Indian Foot, which are the only feet supported as units by my current version of PROJ, hence the explicit scaling factor in the PROJ string.

Note that if you use sf spatial objects instead of sp, then once you've set the CRS it will compute areas and give the units of area, but make sure you are using the latest sf package since Edzer only fixed a bug that ignores +to_meters parameters today (after I reported it because I was trying sf with your projection).

Related Question