[R-sig-Geo] misalignment between data projected in R and ArcGIS (Helen Sofaer)

Ralf Schäfer senator at ecotoxicology.de
Thu Sep 26 18:59:14 CEST 2013


Dear Helen, dear all,

it seems that indeed R is the source of the error.

For the points in Decimal Degrees
> -108.58333   48.75000


I checked with GDAL tools in Shell:

gdaltransform -s_srs EPSG:4326 -t_srs '+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0'
-108.58333 48.75000

and obtain
-871738.467677399 1091090.13492387

System information:
gdaltransform --version
GDAL 1.10.0, released 2013/04/24
uname
Darwin
(OS X 10.8.5)


whereas in R I obtain with:
library(sp)
dat <- SpatialPoints(matrix(data = c(-108.58333, 48.75000), ncol=2))
dat at proj4string  <- CRS("+proj=longlat +datum=WGS84")
library(rgdal)
spTransform(dat, CRS("+proj=aea +lat_1=20 +lat2=60 +lat0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))

 -1207845   5150590
as Helen already has pointed out.

So something seems to go wrong in spTransform - but others are probably better in checking this - or whether we are somehow confusing sth.

I post my system information 
R version 3.0.2 (2013-09-25)
Platform: x86_64-apple-darwin10.8.0 (64-bit)
locale:
[1] de_DE.UTF-8/de_DE.UTF-8/de_DE.UTF-8/C/de_DE.UTF-8/de_DE.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] rgdal_0.8-11 sp_1.0-13   

loaded via a namespace (and not attached):
[1] grid_3.0.2      lattice_0.20-23 tools_3.0.2 

and rgdal gives:
rgdal: version: 0.8-11, (SVN revision 479M)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08
Path to GDAL shared files: /Users/ralfs/Library/R/3.0/library/rgdal/gdal
Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
Path to PROJ.4 shared files: /Users/ralfs/Library/R/3.0/library/rgdal/proj

Best regards
Ralf


Am 26.09.2013 um 18:05 schrieb Helen Sofaer <helen at rams.colostate.edu>:

> Thanks for taking the time to look at it Ralf.
> I wish that was the problem (I had transformed from WGS 84). Here's a
> simpler example that just relies on the future climate data, since I
> realized there's no reason to bring in that state boundary data beyond
> to show that both programs were consistent with themselves for
> different data.
> What I find is that if I project the data in R I get very different
> results than if I project the data in ArcCatalog. The reason this is
> so surprising to me is that I'm defining my proj4 string to match what
> is coming from GIS.
> 
> I may very well be messing this up somehow within ArcGIS, rather than
> within R; it does seem like I've made a simple but major mistake
> somewhere along the way (and I admit I'm more comfortable in R, which
> is why I'm going through this to begin with). Can anyone else with GIS
> (I have 10.1) replicate this issue? The decimal degrees used in my
> example are at the end of the code, but I'd guess that one would see
> the same thing with any other point locations.
> 
> Thanks again,
> Helen
> 
> 
> # Simple version:
> library(sp)
> library(rgdal)
> 
> # just take 10 rows (see earlier code for importing data, or copy
> long/lat from end)
> set.seed(123)
> CGCM.Sample10 = CGCM.A2.May2047[sample(nrow(CGCM.A2.May2047), 10), ]
> write.csv(CGCM.Sample10, "CGCM.Sample10.csv") # for bringing into ArcGIS
> 
> # convert to SPDF, assuming NAD83
> CGCM.Sample10 = SpatialPointsDataFrame(data = CGCM.Sample10,
>            coords = CGCM.Sample10[, c("long_dd", "lati_dd")],
>            proj4string = CRS(as.character("+proj=longlat +datum=NAD83")))
> # project to Albers
> CGCM.Sample10.Albers = spTransform(CGCM.Sample10, CRS("+proj=aea
> +lat_1=20 +lat2=60 +lat0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83
> +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
> bbox(CGCM.Sample10.Albers)
> #             min       max
> # long_dd -1617766 -136530.6
> # lati_dd  4702653 5158505.5
> 
> # Import same points from ArcGIS; CGCM.Sample10.csv was brought into
> ArcCatalog as XY table in NAD83 and converted to shapefile (i.e.
> without projected to Albers within GIS):
> CGCM.Sample10.viaGIS.NAD83 =
> readOGR("C:/Helen_PPR_GIS/ForPondModel/SpatialExtent",
> "CGCM_Sample10_viaGIS_NAD83")
> print(proj4string(CGCM.Sample10.viaGIS.NAD83))
> # [1] "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"
> 
> # if I project these data in R, I get the same results as before:
> CGCM.Sample10.Albers.NotProjectedinGIS =
> spTransform(CGCM.Sample10.viaGIS.NAD83, CRS("+proj=aea +lat_1=20
> +lat2=60 +lat0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m
> +no_defs +ellps=GRS80 +towgs84=0,0,0"))
> bbox(CGCM.Sample10.Albers.NotProjectedinGIS)
> #               min       max
> # coords.x1 -1617766 -136530.6
> # coords.x2  4702653 5158505.5
> 
> ## The problem arises when I import data that have already been
> projected within ArcGIS
> # Using 'North America Albers Equal Area Conic' projection in GIS
> CGCM.Sample10.viaGIS.Albers =
> readOGR("C:/Helen_PPR_GIS/ForPondModel/SpatialExtent",
> "CGCM_Sample10_viaGIS_Albers")
> print(proj4string(CGCM.Sample10.viaGIS.Albers))
> # [1] "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0
> +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
> bbox(CGCM.Sample10.viaGIS.Albers)
> #                 min       max
> # coords.x1 -1238062.2 -100922.4
> # coords.x2   481037.2 1091090.2
> 
> # Full coordinates projected in R:
> cbind(CGCM.Sample10.Albers at coords, CGCM.Sample10.Albers at data$long_dd,
> CGCM.Sample10.Albers at data$lati_dd)
> # Output (Albers X, Y, Decimal degrees X, Y):
>                                            long_dd       lati_dd
> [1,] -1207845.6   5150590     -108.58333   48.75000
> [2,]  -386539.6    4920224     -100.00000   46.33333
> [3,] -1010183.5   5072255     -106.50000   47.91667
> [4,]  -231920.7    5142912     -98.41667    48.91667
> [5,]  -136530.6    5006086     -97.41667    47.33333
> [6,] -1617766.1   4801139     -112.66667   44.58333
> [7,]  -818944.6    4702653     -104.41667  43.83333
> [8,]  -216357.0    5078695     -98.25000    48.16667
> [9,]  -775495.3    5158506     -104.08333  49.00000
> [10,]  -934670.2   4909220     -105.66667  46.08333
> 
> # Full coordinates projected in GIS:
> cbind(CGCM.Sample10.viaGIS.Albers at coords,
> CGCM.Sample10.viaGIS.Albers at data$long_dd,
> CGCM.Sample10.viaGIS.Albers at data$lati_dd)
> # Output (Albers X, Y, Decimal degrees X, Y):
> [1,]  -871738.7   1091090.2    -108.58333 48.75000
> [2,]  -289861.4  754322.3     -100.00000 46.33333
> [3,]  -738909.7  975962.3     -106.50000 47.91667
> [4,]  -167394.6 1054895.9     -98.41667 48.91667
> [5,]  -100922.4  866975.6     -97.41667 47.33333
> [6,] -1238062.2  650570.7    -112.66667 44.58333
> [7,]  -635489.0  481037.2     -104.41667 43.83333
> [8,]  -157948.0  966332.5     -98.25000 48.16667
> [9,]  -558453.2 1086390.0    -104.08333 49.00000
> [10,]  -702497.8  754532.6    -105.66667 46.08333
> 
> On Thu, Sep 26, 2013 at 4:34 AM, Ralf Schäfer <senator at ecotoxicology.de> wrote:
>> Dear Helen,
>> 
>> I do not have access to ArcGIS - so I can not even open the file (I can
>> unzip the .lpk but not open the .sdc file). However, the package information
>> says that the data is in WGS 84 - so if you imported into R you need to
>> assign WGS 84 first and then convert and this should work.
>> From memory, ArcGIS or QGIS always assign a CRS on opening so make sure it
>> is the correct one before projection.
>> From my experience, in 99% of cases where students struggled with the CRS it
>> was an issue with QGIS and ArcGIS.
>> 
>> Best regards,
>> 
>> Ralf Schäfer
>> 
>> ------------------------------------------------------------
>> 
>> Prof. Dr. habil. Ralf Bernhard Schäfer
>> Juniorprofessor for Quantitative Landscape Ecology
>> Environmental Scientist (M.Sc.)
>> Institute for Environmental Sciences
>> University Koblenz-Landau
>> Fortstrasse 7
>> 76829 Landau
>> Germany
>> Mail: schaefer-ralf at uni-landau.de
>> Phone: ++49 (0) 6341 280-31536
>> Web: www.landscapecology.uni-landau.de
>> 
>> 
>> 
>> 
>> 
>> 
>> 
>> 
>> 
>> 
>> 
> 
> 
> 
> -- 
> Helen Sofaer
> Postdoctoral Fellow
> Fish Wildlife and Conservation Biology
> Colorado State University



More information about the R-sig-Geo mailing list