[R-sig-Geo] Error when projecting polygons into Azimuthal Equidistant Projection

[bart]

Hi All

First a short introduction to what i'm trying to do, I'm trying to find 
the length of a line going through a point intersecting with the coast 
lines. My thought was to find the right polygon and then re-project it 
to an Azimuthal Equidistant Projection on which an intersection with the 
line could work using gIntersection from the rgeos package.

Some how re-projecting using spTransfrom invalidates the polygon. The 
polygon starts intersecting with itself. Does anyone have any suggestions?

A small example to reproduce the problem:

require(maps)
require(maptools)
require(rgeos)
require(rgdal)
a<-getRgshhsMap(xlim=c(-170, 180),ylim=c( -60,90))
a<-a[1]# get only eurasia
gIsValid(a)
plot(a)
transformed<-spTransform(a, CRS("+proj=aeqd +lon_0=-9 +lat_0=39" ))
gIsValid(transformed)
x<-9.35151 *10^6
y<-3.60185 *10^6
plot(transformed); axis(1); axis(2); points(x,y, col="red")
plot(transformed, xlim=x+c(-10^5, 10^5), ylim=y+c(-10^5, 10^5)); 
axis(1); axis(2); points(x,y, col="red")


I tried interpolating the segments along a great circle route but that 
did not help:

require(geosphere)
ps<-as.data.frame(SpatialPolygons2PolySet(a))
ps<-ps[ps$SID==1,] 
psnew<-as.data.frame(gcIntermediate(as.matrix(ps[-nrow(ps),c("X","Y")]),as.matrix(ps[-1,c("X","Y")]), 
sepNA=TRUE, addStartEnd=TRUE, n=100))
        psnew<-psnew[!is.na(psnew$lat),]
        psnew<-psnew[!duplicated(psnew),]
        names(psnew)<-c("X","Y")
        psnew$PID<-1
        psnew$POS<-1:nrow(psnew)
        spnew<-PolySet2SpatialPolygons(as.PolySet(psnew,projection="LL"))
        gIsValid(spnew)
        sptra<-spTransform(spnew, CRS("+proj=aeqd +lon_0=-9 +lat_0=39" ))
        x<-9.38019 *10^6
        y<-3.55167 *10^6
        plot(sptra); axis(1); axis(2); points(x,y, col="red")
        plot(sptra, xlim=x+c(-10^5, 10^5), ylim=y+c(-10^5, 10^5)); 
axis(1); axis(2); points(x,y, col="red")
        plot(sptra, xlim=x+c(-10^4, 10^4), ylim=y+c(-10^4, 10^4)); 
axis(1); axis(2); points(x,y, col="red")
        gIsValid(sptra)

Thanks in advance for any tips

Bart

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo