[R-sig-Geo] misalignment between data projected in R and ArcGIS (Helen Sofaer)
Roger Bivand
Roger.Bivand at nhh.no
Thu Sep 26 20:25:49 CEST 2013
Summary: typo not error in OSX binary. Detail below:
On Thu, 26 Sep 2013, Roger Bivand wrote:
> 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
Sorry, no, this is a simple typo, +lat0=40 for +lat_0=40. With the latter,
OK, no issue with PROJ versions. If you set +lat_0=0 - what the typo gives
- you get the erroneous value. Please check on OSX (Helen is on Windows, I
think, from a file path).
> aea1 <- "+proj=aea +lat_1=20 +lat2=60 +lat_0=0 +lon_0=-96 +x_0=0 +y_0=0
+datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
> project(pt, aea1)
[,1] [,2]
[1,] -1207845 5150590
>> +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.
There is no reason to suspect the CRAN OSX binary of error!
>
> 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.
There is no reason to suspect the CRAN OSX binary of error!
>
> 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.
This may make sense anyway ...
Roger
>
> 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
^^^^^^^^ ^^^^^^^^ more issues
>>> +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