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

[bart]

Thank you very much for your answer Roger, but the problem is that i 
want to do this operation for a regular grid around the world. Meaning i 
will encounter the error in some places. As far as i understand the 
projection the distance and direction to other locations is only valid 
if lat_0 and lon_0 are the location of interest. So avoiding this 
projection is not an option.

Bart
On 08/15/2011 06:51 PM, Roger Bivand wrote:
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@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo



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