R help - Sep 2012 - why variations in accuracy between R to ARCGIS for the same point reprojection?

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/

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:
> 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/
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595