# Understanding the unit of area for standard deviational ellipse

I am using the `aspace` package to calculate the standard deviational ellipse of a set of points. Test shapefile with 5 points zipped here.

``````library(maptools)
library(sp)
library(aspace)
library(rgdal)
coords <- shpfile@coords
calc_sde(points=coords)
plot_sde(plotnew = TRUE)
r.SDE\$Area.sde
``````

The area is 2824712729.

I think I am getting tripped up by the coordinate system:

``````coords
#  coords.x1 coords.x2
#0   3667152   8861968
#1   3713150   8878606
#2   3714374   8843128
#3   3683545   8834809
#4   3690396   8854383
``````

The shapefile is set to `+proj=longlat +datum=WGS84 +no_defs`.

1. What are the units for this area?
2. How can I project the coordinates in decimal degrees?

• Those coordinates are not lat/long decimal degrees. You need to find out the actual coordinate system. I would imagine that the projection units are in meters. Mar 19, 2017 at 3:30
• In your picture your units are meters, so your units for area would be square meters. `2824712729m^2` or `2824.7km^2`. Looking at your grid that area looks about right.
– Midavalo
Mar 19, 2017 at 5:14
• Your title has "SDE" ((which hasn't existed in a long time) and your question body references "SHP" (which is a suffix, not a format) -- these are not the same thing. Mar 19, 2017 at 11:25
• @JeffreyEvans I think it's WGS84. Is there a way to confirm? Mar 19, 2017 at 11:42
• @Kazuhito this is what I needed for Question 2. As JeffreyEvans and Midavalo noted, the units are meters. Still looking into whuber's note about the incorrect ellipse. Mar 19, 2017 at 19:41

``````library(maptools); library(sp); library(aspace); library(rgdal)
proj4string(shpfile) <- CRS("+init=epsg:32736")  # UTM36S
coords <- shpfile@coords
calc_sde(points=coords)
plot_sde(plotnew = TRUE)
r.SDE\$Area.sde
# [1] 2824692311 (m2)

shpfile_lonlat <- spTransform(shpfile, CRS("+proj=longlat +ellps=WGS84"))
LonLat <- shpfile_lonlat@coords
calc_sde(points=LonLat)

# export SDE as shapefile
shp <- shapefiles::convert.to.shapefile(sdeloc,sdeatt,"id",5)
write.shapefile(shp, "./SDE_shape", arcgis=T)

# read SDE shapefile
proj4string = CRS("+proj=longlat +ellps=WGS84"))
plot(Ellips, border= "blue", axes= TRUE)
points(shpfile_lonlat)

# calculate area on WGS84 ellipsoid
library(geosphere)
areaPolygon(Ellips, a=6378137, f=1/298.257223563) # WGS84
# [1] 2216829475 (m2)
``````

As `r.SDE\$Area.sde` seems to calculate planimetric area, I tried ellipsoidal area on WGS84 using {geosphere} package.

• this is great, @Kazuhito Mar 22, 2017 at 1:51
• Hi @EricGreen Correction to the code with apologies. In `write.shapefile()` function `"SDE_shape.shp"` should have been `"SDE_shape"` without .shp (edited the answer). Apparently you have fixed the code yourself :) Mar 23, 2017 at 1:46
• Yes, thanks. I found I needed to make that small change. Mar 23, 2017 at 1:49