Skip to content

Error when projecting polygons into Azimuthal Equidistant Projection

4 messages · Kranstauber, Bart, Roger Bivand

Original
Kranstauber, Bart
# 
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?

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
Roger Bivand
#
On Mon, 15 Aug 2011, bart wrote:

            

      
      
      
    

        
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

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Kranstauber, 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:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Tue, 16 Aug 2011, bart wrote:

        

            

      
      
      
    

        
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