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

Thu Sep 26 18:05:50 CEST 2013

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