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

Roger Bivand Roger.Bivand at nhh.no
Thu Sep 26 20:01:25 CEST 2013


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


More information about the R-sig-Geo mailing list