Skip to content

why variations in accuracy between R to ARCGIS for the same point reprojection?

3 messages · Brian Ripley, Camilo Mora

Camilo Mora
# 
Hi everyone,

I wonder if anyone knows the reason why the outputs of the same  
reprojection in r and arcgis are different?. The magnitude of the  
change can be up to 40 km in the poles.
Basically, I have a database of points equally separated by one degree  
over the globe.
In ARCGIS,  I am projecting the data in GCS-WGS-1984 and then  
reprojected it to Berhmann to ensure equal area distribution of the  
points.
In R, I am using:
spPoint <-  
SpatialPoints(coords=coordinates(Data),proj4string=CRS("+proj=longlat  
+datum=WGS84"))
and then reprojecting it to Berhmann with:
spPointReprj=spTransform(Data,CRS("+proj=cea +lon_0=0 +lat_ts=30  
+x_0=0 +y_0=0 +ellps=WGS84"))

If I put the two outputs of the reprojections in the same map, they  
are off by few meters in the tropics by up to 40km in the poles.

I decided to investigate how the reprojections are done and my  
calculations are different from both R and ArcGis:

First, I calculated the radious of the Earth in two different ways:
=Re * (1 - e^2)/ (1 - e^2 *SIN(RADIANS(Latitude))^2)^(3/2)
=Re * SQRT(1 - e^2) / (1 - e^2 * (SIN(RADIANS(Latitude)))^2)

where Re is the radius of the Earth at the tropics(6378km) and e is  
the eccentricity of the ellipsoid (0.081082).

According to several books I used, the position of a point in the  
Y-axis in the Berhmann projection could be estimated as:
=Re*(SIN(RADIANS(Latitude))/COS(RADIANS(Spll)))
where Spll is the standard parallel, which in the Berhmann's projection is 30.
Unfortunately, the resulting output, although similar in shape to the  
outputs in R and Arcgis, is still not quite the same. Any thoughts why  
these differences in supposedly the same calculations?

Any input will be greatly appreciated,

Thanks,

Camilo




Camilo Mora, Ph.D.
Department of Geography, University of Hawaii
Currently available in Colombia
Phone:   Country code: 57
          Provider code: 313
          Phone 776 2282
          From the USA or Canada you have to dial 011 57 313 776 2282
http://www.soc.hawaii.edu/mora/
Brian Ripley
# 
There is no 'reprojection' in R (which is upper case).  Please attribute 
blame correctly.

You seem to be talking about some contributed addon package, not 
specified.  But I think you should be asking this on R-sig-geo.
On 02/09/2012 19:24, Camilo Mora wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Camilo Mora
    

    

      
      #
    

  


  

      
I apologize, I did not intend to blame anyone. I actually thought  
there may be some underlying differences in the formulas used. The  
patterns overall are very similar but they are not that precise to one  
another. The function I use in r is called spTransform from the  
package "rgdal".

Again, my guess is that it has to do something with the functions  
used. As an example, I found two different ways to measure the radius  
of the Earth. So perhaps there is an additional parameter that  
differentiates the scripts from R and ArcGis?.

Sorry again for the confusion and thanks,

Camilo
p.s. I always pride R. It is just hard to image academic life without it...

Camilo Mora, Ph.D.
Department of Geography, University of Hawaii
Currently available in Colombia
Phone:   Country code: 57
          Provider code: 313
          Phone 776 2282
          From the USA or Canada you have to dial 011 57 313 776 2282
http://www.soc.hawaii.edu/mora/



Quoting Prof Brian Ripley <ripley at stats.ox.ac.uk>: