Dear list,

I've been having trouble trying to generate a polygon from points crossing the international date line (or antemeridian) and I was wondering if anyone has come up with a solution to the problem?

Attached are the locations from a tagged seabird foraging  off the NE coast of New Zealand, my code is as follows...

Any help is much appreciated.

Kind regards,

Ross

Dr Ross Dwyer | School of Biological Sciences | University of Queensland | Email: ross.dwyer@uq.edu.au<mailto:ross.dwyer@uq.edu.au>


###############################

library(maps)
library(adehabitatHR)
library(plotKML)
library(shapefiles)
library(rgeos)

#Create the data
Latitude <- c(-37.15,-38.04,-37.92,-36.78,-36.6,-35.50,-34.67,-35.63,-36.88)
Longitude <- c (-160.02,-158.58, -172.39, -169.95, -160.01,  178.69, 169.88, -170.25, -172.81)
inputfile <- data.frame(Latitude, Longitude)

mylocs<- inputfile
coordinates(mylocs) <- ~Longitude+Latitude
proj4string(mylocs) <- CRS("+proj=longlat +datum=WGS84")

#Plot the data
map('world')
plot(mylocs,add=T,col=2)
mymcp <- mcp(mylocs,percent=100)
plot(mymcp,add=T)
# Or in Google Earth
kml(mymcp,file='mynewmcp.kml')

#You'll see that the polygon is stretched across the globe.

#What I want to happen is...
#step 1 add 360 to negatives to make all xys on same plane
inputfile$newLong <- ifelse(inputfile$Longitude < 0,inputfile$Longitude+360,inputfile$Longitude)

mylocs2<- SpatialPoints(inputfile[,c("newLong","Latitude")])
map('world2', xlim = c(130, 300))
plot(mylocs2,add=TRUE,col=2)
mymcp2 <- mcp(mylocs2,percent=100)
plot(mymcp2,add=T)

# Ultimately what I would like to do is extract the area of this polygon (after reprojecting into the relevant coordinate system) and to save the polygon object as a .kml for visualisation in Google Earth.
#  I've tried splitting the polygon across the antemeridian  but run into difficulties when there are multiple polygons within a SPDF object.

# Find where lines intersect mcp polygon
L1 = Line(cbind(c(180,180),c(-90,90)))
Ls1 = Lines(list(L1), ID = "am")
antemeridlines = SpatialLines(list(Ls1))
plot(antemeridlines,add=T,col=3)
myintersect <- gIntersection(antemeridlines,mymcp2)
myintersectY <- coordinates(myintersect)[[1]][[1]][,2]
myintersectY

###############################



	[[alternative HTML version deleted]]

