2

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)
shpfile <- readShapePoints("testSHP.shp")
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?

enter image description here

Improve this question
11
  • 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.
    – Jeffrey Evans
    Commented 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 ♦
    Commented 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
    Commented Mar 19, 2017 at 11:25
  • @JeffreyEvans I think it's WGS84. Is there a way to confirm?
    – Eric Green
    Commented 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
    Commented Mar 19, 2017 at 19:41

1 Answer 1

Reset to default
4 
library(maptools); library(sp); library(aspace); library(rgdal)
shpfile <- readShapePoints("testSHP.shp")
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
Ellips <- readShapeSpatial("./SDE_shape.shp", 
                            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.

Improve this answer
3
  • this is great, @Kazuhito
    – Eric Green
    Commented 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
    Commented Mar 23, 2017 at 1:46
  • Yes, thanks. I found I needed to make that small change.
    – Eric Green
    Commented Mar 23, 2017 at 1:49

Your Answer

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

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