misalignment between data projected in R and ArcGIS (Helen Sofaer)
On Thu, 26 Sep 2013, Ralf Sch?fer wrote:
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.
The missing information is about the way in which rgdal was installed. For my rgdal, installed from source and dynamically linked to the PROJ library (GDAL also links to the PROJ library, so PROJ is the crucial component), I see: $ proj +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 -871738.47 1091090.13 So OK, and:
library(rgdal)
Loading required package: sp rgdal: version: 0.8-11, (SVN revision 479M) Geospatial Data Abstraction Library extensions to R successfully loaded Loaded GDAL runtime: GDAL 1.10.1, released 2013/08/26 Path to GDAL shared files: /usr/local/share/gdal Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480] Path to PROJ.4 shared files: (autodetected)
pt <- matrix(c(-108.58333, 48.75000), nrow=1) project(pt, "+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")
[,1] [,2]
[1,] -871738.5 1091090
a <- project(pt, "+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")
b <- spTransform(SpatialPoints(pt, proj4string=CRS("+init=epsg:4326")),
CRS("+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"))
all.equal(a, coordinates(b), check.attributes=FALSE)
[1] TRUE the same in R. The failing component is then the underlying PROJ library. Ralf is using OSX, and if both Ralf and Helen are using the CRAN OSX rgdal binary, this may be the explanation, but of course will not be if rgdal and PROJ were installed in a different way. I have checked that the CRAN Windows rgdal binary works correctly, using PROJ 4.7.1, linked statically, on 64-bit and 32-bit. Please describe how you installed rgdal, and if it turns out to be the OSX CRAN binary, we'll address that - if however you installed PROJ in some other way then rgdal from source, we'll need to check those routes instead. This also suggests that (part of) the PROJ.4 test suite should be added to rgdal to check the integrity of the projection and datum transformation functions. Thanks for clear examples, Roger
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
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand Department of Economics, NHH Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand at nhh.no