On Tue, 16 Aug 2011, bart wrote:
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.
The underlying projection uses PROJ.4, as I'm sure you are aware. Consequently, you should probably choose a different projection, since +proj=aeqd will always fail in this way for your choice of +lon_0.
Could you extract the intersection points in rgeos in geographical coordinates, and, after finding the appropriate pairs, measure the distance with spDistsN1(..., longlat=TRUE)?
Roger
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. RogerA 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_______________________________________________ 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
