[R-sig-Geo] rbind of SpatialPolygonDataFrame objects is different with or without rgdal package

[Roger Bivand]

On Mon, 10 Oct 2011, MacQueen, Don wrote:

> It appears that when I rbind() a pair of SpatialPolygonDataFrame objects
> when the rgdal package has been loaded the resulting projection string
> gets an embedded space character, whereas without rgdal it does not. The
> space character breaks further use of rbind(). I don't think it matters
> which one is first in the search path (I've tried it both ways).
>

This is a result of our choice a long time ago not to package PROJ.4 with 
sp but with rgdal (GDAL depends on PROJ.4, so it had to be there anyway). 
Consequently, code checking the validity of the PROJ.4 string used to set 
the CRS object only does a cursory job when only sp is loaded, but calls 
out to PROJ.4 when rgdal is loaded. PROJ.4 replies and a space gets 
prepended (which might be fixable), and in your case also appends a 
+towgs84=0,0,0, which isn't fixable. If rgdal is loaded and the CRS object 
is validated, it may be normalised in this way by PROJ.4 using 
pj_get_def() in that library, which appears to insert the initial space.

rbind() in sp checks the projection, so makes the change when and only 
when rgdal is loaded. spRbind() in maptools - arguably wrongly - does not 
check the projection, so works the same on pairs but not dotlists of 
objects whether or not rgdal is loaded.

So the workaround may be to load rgdal and normalise all the affected 
objects with CRSargs().

It isn't obvious how to do this right, I'm afraid.

Hope this clarifies,

Roger

>
>
> Here's an example.
>
> I have three simple spatial polygon data frame objects; each consists of a
> single polygon. They have the same projection. I can send them as an
> attachment if requested.
>
> ## in a new R session (objects created in an earlier session):
>
>> require(sp)
> Loading required package: sp
>
>> proj4string(ex1.sp) == proj4string(ex2.sp)
> [1] TRUE
>
>> proj4string(ex1.sp) == proj4string(ex3.sp)
> [1] TRUE
>
>> tmp1 <- rbind(ex1.sp, spChFIDs(ex2.sp,'ex2'))
>
>> proj4string(ex1.sp) == proj4string(tmp1)
> [1] TRUE
>
>> sessionInfo()
> R version 2.13.1 (2011-07-08)
> Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
>
> locale:
> [1] C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] sp_0.9-88
>
> loaded via a namespace (and not attached):
> [1] grid_2.13.1     lattice_0.19-30
>
>
>> require(rgdal)
> Loading required package: rgdal
> Geospatial Data Abstraction Library extensions to R successfully loaded
> Loaded GDAL runtime: GDAL 1.8.0, released 2011/01/12
> Path to GDAL shared files:
> /Library/Frameworks/R.framework/Versions/2.13/Resources/library/rgdal/gdal
> Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009, [PJ_VERSION: 470]
> Path to PROJ.4 shared files:
> /Library/Frameworks/R.framework/Versions/2.13/Resources/library/rgdal/proj
>
>
>
>> tmp2 <- rbind(ex1.sp, spChFIDs(ex2.sp,'ex2'))
>
>> proj4string(ex1.sp) == proj4string(tmp2)
> [1] FALSE
>
>> proj4string(ex1.sp)
> [1] "+proj=lcc +lat_1=38.43333333333333 +lat_2=37.06666666666667
> +lat_0=36.5 +lon_0=-120.5 +x_0=2000000.0001016 +y_0=500000.0001016001
> +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs"
>
> ## note the leading space character
>> proj4string(tmp2)
> [1] " +proj=lcc +lat_1=38.43333333333333 +lat_2=37.06666666666667
> +lat_0=36.5 +lon_0=-120.5 +x_0=2000000.0001016 +y_0=500000.0001016001
> +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs
> +towgs84=0,0,0"
>
>> tmp3 <- rbind(tmp2, spChFIDs(ex3.sp,'ex3'))
> Error in checkCRSequal(dots) : coordinate reference systems differ
>
>
>> sessionInfo()
> R version 2.13.1 (2011-07-08)
> Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
>
> locale:
> [1] C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] rgdal_0.7-1 sp_0.9-88
>
> loaded via a namespace (and not attached):
> [1] grid_2.13.1     lattice_0.19-30
>>
>
>
>
>> str(ex1.sp)
> Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
>  ..@ data       :'data.frame': 1 obs. of  1 variable:
>  .. ..$ name: Factor w/ 1 level "example1": 1
>  ..@ polygons   :List of 1
>  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
>  .. .. .. ..@ Polygons :List of 1
>  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
>  .. .. .. .. .. .. ..@ labpt  : num [1:2] 6223669 2081123
>  .. .. .. .. .. .. ..@ area   : num 20287432
>  .. .. .. .. .. .. ..@ hole   : logi FALSE
>  .. .. .. .. .. .. ..@ ringDir: int 1
>  .. .. .. .. .. .. ..@ coords : num [1:8, 1:2] 6221449 6226373 6225988
> 6223664 6223594 ...
>  .. .. .. ..@ plotOrder: int 1
>  .. .. .. ..@ labpt    : num [1:2] 6223669 2081123
>  .. .. .. ..@ ID       : chr "1"
>  .. .. .. ..@ area     : num 20287432
>  ..@ plotOrder  : int 1
>  ..@ bbox       : num [1:2, 1:2] 6220680 2078695 6226373 2083511
>  .. ..- attr(*, "dimnames")=List of 2
>  .. .. ..$ : chr [1:2] "x" "y"
>  .. .. ..$ : chr [1:2] "min" "max"
>  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
>  .. .. ..@ projargs: chr "+proj=lcc +lat_1=38.43333333333333
> +lat_2=37.06666666666667 +lat_0=36.5 +lon_0=-120.5 +x_0=2000000.0001016
> +y_0=500000.00010160"| __truncated__
>
>
>
>
>
>

-- 
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