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

[Roger Bivand]

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.
>> 
>> 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
>>> 
>> 
>
> _______________________________________________
> 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