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