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

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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