gArea() units

I would strongly advise against using EPSG:32640 for any region that
is not nearby to longitude 57E (where "nearby" means within about 6
degrees of longitude (but I'm very wary of providing advice here it's
something I'm not expert in).

Projection choice is really important, and UTM is simply not
appropriate for many applications (in the same way that no projection
is appropriate for all applications). It's certainly not appropriate
to use a UTM zone that is far from the region of study.

Look at the numbers that come out:

## my best-guess
gArea(spTransform(myPolygonSpatial, CRS("+proj=laea +lon_0=-1.92 +
+lat_0=55.7 +datum=WGS84")))/1e6
[1] 5.796854

## UTM that is for a region 56 zones away(!)
gArea(spTransform(myPolygonSpatial, CRS("+init=epsg:32640")))/1e6
[1] 7.475327

A better choice if you must use UTM, would be the right zone:

gArea(spTransform(myPolygonSpatial, CRS("+init=epsg:32631")))/1e6
[1] 5.805833

As Roger points out there is areaPolygon() in geosphere, a good way to
do a check on whether we are in the right ball-park:

areaPolygon(myPolygonSpatial)/1e6
[1] 5.782704

Another thing to worry about if you are reprojecting to get the
calculation properties that you need is the amount of "segmentation"
in the geometry. A "straight-line" in one projection defined by only
two end points can be a rather curved shape in another projection, and
so the shape itself can only be preserved by inserting extra
"redundant" points first. That's not a problem in this particular
example for such a small area on the Earth, though the wrong UTM zone
really is a problem.  (Sometimes itt just depends on the specifics to
an example. . . )


I don't want to get into advising about projections without doing a
lot more work of my own to ensure my advice is sound, but please be
aware that this is something that needs quite a lot of thought by the
user, or alternatively---expert local advice.


Cheers, Mike.

Hi Ross,

you result is in square degrees, you will need to project your data to
planar coordinates.

You could do for example:

library(rgeos)
library(rgdal)
library(sp)

myPoly <- spTransform(myPolygonSpatial, CRS("+init=epsg:32640"))

# the area in m2
 gArea(myPoly)

# or in km2
gArea(myPoly) / 1e6

HTH

Johannes




On Tue, May 14, 2013 at 4:08 PM, Ross Ahmed <rossahmed at googlemail.com>wrote:

Hi all

The documentation for gArea() says the units of area are based on the
current projection. The polygon below is in longlat, and in my primative
understanding, the units should be square metres. The output is
 0.000828375
but I know the approximate area of this polygon is 4.45 square kilometres.
Using gArea(), how can I calculate the area of the polygon below in units
of
square kilometres?

library(sp); library(rgeos)
polygonCoords <- matrix(c(-1.9450, -1.9075, -1.9075, -1.9450, -1.9450,
55.72476, 55.72476, 55.70267, 55.70267, 55.72476), ncol=2)
p = Polygon(polygonCoords)
myPolygon = Polygons(list(p),1)
myPolygonSpatial = SpatialPolygons(list(myPolygon))

proj4string(myPolygonSpatial) <- CRS("+proj=longlat +datum=WGS84
+ellps=WGS84 +towgs84=0,0,0")

gArea(myPolygonSpatial)
# output is [1] 0.000828375
# area in square kilometres is circa 4.45 square kilometres


Thanks
Ross





        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

--

## ## Contact:
## Johannes Signer, MSc
## PhD-Candidate -- Wildlife Management
## Dept. of Forest Zoology & Forest Conservation
## University of Goettingen
## Buesgenweg 3 (room 85)
## 37077 Goettingen, Germany
## e: jmsigner at gmail.com
## p: +49 (0)551 39 3583
## m: +49 (0)176 61 962856
## skype: j.m.signer
## ##

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Michael Sumner
Hobart, Australia
e-mail: mdsumner at gmail.com