0

我在 R 计算动物的活动范围时遇到了很多麻烦。我想一旦我制作了一个家庭范围（如果我做得正确的话）计算面积会很容易，但没有

我已经粘贴了一些我一直在尝试的代码。我想知道有人有任何见解吗？

# Load package

    library(adehabitat)

#Load file Frodo

    dd <- read.csv(file.choose(), header = T)

# Plot the home range

    xy <- dd[,c("X","Y")]
    id <- dd[,"name"]
    hr<- mcp(xy,id,percent=95)
    plot(hr)
    points(xy[xy$id=="frodo",])

#Great. Home range produced.  Now calculate area

    area <- mcp.area(xy, id,percent = 95),

# Result 2.287789e-09 Ha.  Way to small.  Maybe it doesnt like Lat / Long.  

# Will try and convert coordinates  into M or Km

# Load map project

    library(mapproj)

    x<-mapproject(t$X,t$Y,projection="mercator")

# Its converted it to something but its not M's  or Km's.  
# I'll try and run it anyway

    xy <- x[,c("X","Y")]

# incorrect number of dimensions

# Ill try Project 4
library(proj4)

    xy <- dd[,c("X","Y")]
    tr <- ptransform(xy/180*pi, '+proj=latlong +ellps=sphere',
                 '+proj=merc +ellps=sphere')
    View(tr)

# There seems to be a Z column filled with 0's. 
# It that going to affect anything? 
# Let's look at the data
    plot(tr)

# Looks good, Lets try and create a home range

    xy <- tr[,c("x","y")]
#  'incorrect number of dimensions'

不知道是什么问题。不知道我是在正确的轨道上还是在做完全错误的事情

4

1 回答 1

2

为了计算面积，您需要投影坐标系中的点（长/纬度的面积只是度数单位）。您使用的投影类型将对生成的区域产生很大影响。例如，墨卡托投影扭曲了远离赤道的区域——您可能想要查看您所在位置的最佳等面积投影。我将回答您问题的编程部分，一旦您找到要使用的正确预测，您就可以将它们替换掉。

require(sp)
require(rgdal)
orig.points <- dd[,c("X","Y")]

# geographic coordinate system of your points
c1 <- CRS("+proj=latlong +ellps=sphere") 

# define as SpatialPoints
p1 <- SpatialPoints(orig.points, proj4string=c1) 

# define projected coordinate system of your choice, I am using the one you 
# defined above, but see:     
# http://www.remotesensing.org/geotiff/proj_list/mercator_1sp.html
# to make sure your definition of the mercator projection is appropriate
c2 <- CRS("+proj=merc ellps=sphere")

p2 <- spTransform(p1, c2) # project points
# convert to Polygon (this automatically computes the area as an attribute)
poly <- Polygon(p2)
poly@area #will print out the area
于 2013-08-29T16:11:56.977 回答