# Getting area of shapefile polygons in square kilometers using R?

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)
 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.

• crs(x) <- CRS("+init=epsg:2314") does not change the projection, it assigns it - and area() will assume longlat if if is NA, was there not a warning to this effect? I don't think this is answerable given the information provided, can you share the shapefile? Oct 7 '19 at 15:55
• I did not get any warning, and unfortunately I cannot share the shapefile for sake of confidentiality. Thank you so much for your assistance though. Oct 7 '19 at 16:06
• Can you show us something about your data? What do you see when you `print(x)`? What's the extent? That will tell us if you have coordinates in lat-long or a projected system in metres. As it is we've got nothing. Oct 7 '19 at 17:38

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
 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).

• Oh wow! I would have never noticed that! Thanks so much!!! Oct 8 '19 at 2:04
• Boom, nice one. Oct 8 '19 at 2:15