[R-sig-Geo] misalignment between data projected in R and ArcGIS (Helen Sofaer)
Ralf Schäfer
senator at ecotoxicology.de
Fri Sep 27 00:40:15 CEST 2013
Thanks, heaps Roger!!!
Have no idea how this happened - as I obviously just copy and pasted into the shell window and into R...
Anyway, now works as suggested!
I obtain the same results in R and rgdal
-871738.5 1091090
So, any difficulties seem due to what is done in ArcGIS...
Cheers
Ralf
Am 26.09.2013 um 20:25 schrieb Roger Bivand <Roger.Bivand at nhh.no>:
> 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