# Why does projection not affect answer calculated by rgeos::gArea?

I am trying to calculate the area of various polygons from a shapefile containing areas of land with mining rights is South Africa.

Through this site and posts like:

I understand that I must project to an equal area projection.

So this is some of my code:

Mineralrights <-readOGR(dsn = "C:/blabla/notrelevant/Shapefiles",layer = "Mineral Rights")
spTransform(Mineralrights, CRS("+proj=aea +lat_1=24 +lat_2=-33 +lat_0=0 +lon_0=24 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"))
gArea(Mineralrights[Mineralrights\$GID == "1710350",])
[1] 0.0004728565
Warning message:
In RGEOSMiscFunc(spgeom, byid, "rgeos_area") :
Spatial object is not projected; GEOS expects planar coordinates

And even if I don't perform the spTransform, I get the same warning with the same answer of 0.0004728565

While if I use areaPolygon from the geosphere package, I get a more believable answer:

areaPolygon(Mineralrights[Mineralrights\$GID == "1710350",])
[1] 5130055

Which corresponds much better to two other methods I used to calculate the area (With a polygon measurement tool in Google Earth Pro and a freeware GIS tool called Mapwindow 5) yielding 5,14km^2 and 5.136km^2 respectively.

Why does the answer from gArea stay the same even if I transform to an equal area projection which I think should yield an answer in m^2? (due to the +units=m in the CRS)

• You are executing but not saving the spTransform step. The information you don't share, but is relevant (despite your boilerplate) is what is the CRS of the original read. You want Mineralrights <- spTransform(Mineralrights, CRS("+proj=aea +lat_1=24 +lat_2=-33 +lat_0=0 +lon_0=24 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs")) then run gArea() on that. geosphere is automatically converting from longlat to something-Cartesian-in-metres. Without transforming, gArea is giving Cartesian in longlat. geosphere will in general be more accurate, since not all CRS are area preserving everywhere. – mdsumner Dec 7 '17 at 12:59
• @mdsumner Such a rookie mistake! (facepalm) I just ran it again and got 5129877 for the area, which is exactly what I wanted. For added information, areaPolygon gave 5130055 for the same area. Please add your comment as an answer. – Waldo vdM Dec 7 '17 at 13:22