[R-sig-Geo] misalignment between data projected in R and ArcGIS (Helen Sofaer)
Ralf Schäfer
senator at ecotoxicology.de
Thu Sep 26 18:59:14 CEST 2013
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
More information about the R-sig-Geo
mailing list