3

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.

3
  • 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?
    – mdsumner
    Oct 7, 2019 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.
    – Outlier
    Oct 7, 2019 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.
    – Spacedman
    Oct 7, 2019 at 17:38

1 Answer 1

1

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.