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

shpfile <- readShapePoints("testSHP.shp")
coords <- shpfile@coords
plot_sde(plotnew = TRUE)

The area is 2824712729.

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

#  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?

enter image description here

  • 1
    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.
    – Vince
    Mar 19, 2017 at 11:25
  • @JeffreyEvans I think it's WGS84. Is there a way to confirm?
    – Eric Green
    Mar 19, 2017 at 11:42
  • 1
    @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.
    – Eric Green
    Mar 19, 2017 at 19:41

1 Answer 1

library(maptools); library(sp); library(aspace); library(rgdal)
shpfile <- readShapePoints("testSHP.shp")
proj4string(shpfile) <- CRS("+init=epsg:32736")  # UTM36S
coords <- shpfile@coords
plot_sde(plotnew = TRUE)
# [1] 2824692311 (m2)

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

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

# read SDE shapefile
Ellips <- readShapeSpatial("./SDE_shape.shp", 
                            proj4string = CRS("+proj=longlat +ellps=WGS84"))
plot(Ellips, border= "blue", axes= TRUE)

# calculate area on WGS84 ellipsoid
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
    – Eric Green
    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 :)
    – Kazuhito
    Mar 23, 2017 at 1:46
  • Yes, thanks. I found I needed to make that small change.
    – Eric Green
    Mar 23, 2017 at 1:49

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.