[R-sig-Geo] misalignment between data projected in R and ArcGIS (Helen Sofaer)
Helen Sofaer
helen at rams.colostate.edu
Thu Sep 26 19:55:48 CEST 2013
Thanks Ralf!
This is definitely getting interesting. Someone pointed out to me off
list that it looked like these were decimal degrees in WGS84, rather
than NAD83. (It's not clear from the metadata, and after emailing the
person and not getting a straight answer about the coordinate system,
I had assumed it was the same as their source data, which are in
NAD83.)
When I run it again (all in R) assuming WGS84 for the initial data, I
get the same results as I was getting from GIS (i.e. -871738.7,
1091090.2 for the point Ralf focused on). So from my perspective it
looks like the problem was in my assumption that these were in NAD83.
But, it looks like Ralf didn't go wrong there - his code is already
all for WGS84 - so I'm not sure what's going on.
I had always understood that NAD83 and WGS84 were almost the same, so
I'm surprised that the results I was getting before (and that Ralf is
still getting) are very different when it is done through R, but not
when it is done through GIS.
Helen
On Thu, Sep 26, 2013 at 10:59 AM, Ralf Schäfer <senator at ecotoxicology.de> 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.
>
> 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
>
--
Helen Sofaer
Postdoctoral Fellow
Fish Wildlife and Conservation Biology
Colorado State University
More information about the R-sig-Geo
mailing list