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

[Roger Bivand]

On Mon, 15 Aug 2011, bart wrote:

> 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?

Your choice of +lon_0=-9 is unfortunate and flips Malaysia onto India; 
+lon_0=90 does not lead to meltdown. So it is an artefact of your choice 
of +lon_0 for this projection.

Roger

>
> 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 at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no