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
[email protected]
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: [email protected]
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
