For an interim period, it seems sensible to try to write wrappers for existing packages interfacing spatial data and R, using the new sp classes. This allows the interface packages to continue with their regular codebase, and existing users of these packages to be able to choose when or if they will migrate to sp classes. It also permits us to develop and refine the sp package without dependence on interface packages.
Note that the functionality described here has been included in packages rgdal, naptools, and spgrass6, now all on CRAN. This page is maintained for documentary purposes only, and most of the examples shown here are the examples used in the released packages.
The wrapper packages will be released as R source packages and R Windows binary packages on the R-spatial SourceForge repository (see the download page), and may be checked out for users needing to follow the current CVS code. This document will try to show how to use several of these wrapper packages in conjunction with sp:
Here is an example of using data from a geographical database shipped with the maptools package, here extracting lines rather than polygons. As with polygons, the user may get more in the object returned by the map() function than spefically required. Here, lines beginning within the reqested window are included without clipping to the window, so a small pruning operation may be needed to impose the desired limits.
PBSmapping is a large package on CRAN addressing the mapping needs of fisheries researchers. The package ships with many shoreline data sets, and other data sets for, for example, fisheries management areas. Users may find it helpful to add additional data sets, and the functions in this wrapper are for converting SpatialLines and SpatialPolygons objects from sp to PolySet objects in PBSmapping (thanks to Denis Chabot for this idea).
PolySet
Records : 49004
Contours : 65
Holes : 0
Events : NA
On boundary : NA
Ranges
X : [4.8041934967041, 31.9980335235596]
Y : [58.000804901123, 71.165023803711]
Attributes
Projection : 1
Zone : NULL
Extra columns :
The gpclib package provides facilities for polygon clipping (also available in the PBSmapping package). These may be used for a number of purposes, and this wrapper package will be extended, and may be moved into the sp package as time progresses. SpatialPolygons objects are expected to be constructed from Polygons objects; if an Polygons object contains only one Polygon object, it is expected not to be a hole. If, however, there are multiple Polygon objects in an Polygons object, some may be holes if they are contained within others in the same Polygons object. At the next level down, islands (not holes) may be found, and so on.
Establishing which Polygon members of Polygons objects are holes matters to be able to paint them correctly, because Polygon objects within an Polygons object are painted in decreasing order of gross area. None of the existing functions for importing rings support topology, so polygon clipping is used to determine hole status by taking the intersection of the bounding box of the Polygons object and the component Polygon objects. This sets their hole slots correctly. The checkPolygonsHoles() function is run separately for each Polygons object in an SpatialPolygons object:
Note that both polygon import functions and others, here the asSpatialPolygonsShapes() function, may reorder the Polygons objects; for this reason the IDs should always be checked for coherence, and coherence imposed through matching the data frame row names:
TRUE
TRUE
In this case the read_ShapePoly() imposes the shape ID values (often integers starting at 0) as character IDs for the Polygons objects and the rows of the data frame. Note that these are sorted alphanumerically, giving this output ordering:0 1 2 3 4 5 6 7 8 9
Another popular operation made possible using functions in the gpclib package is aggregation or dissolving of Polygons objects. Here we cut the county centroids (label points) eastings into quartiles, and use these groups to aggregate and dissolve counties into four east-west regions:
[-84.1,-81.2] (-81.2,-79.3] (-79.3,-77.8] (-77.8,-75.9]
When recentering a world map, say to change an "Atlantic" view with longitude range -180 to 180, to a "Pacific" view, with longitude range 0 to 360, polygons crossed by the new offset, here 0/360, need to be clipped into left and right sub.polygons to avoid horizontal scratches across the map. The nowrapSpatialPolygons() function performs this operation using polygon intersection, and the nowrapRecenter() function also recenters the output SpatialPolygons object.
The existing GRASS interface package for R was suited to GRASS 5, and still works with GRASS 5.4. Changes going forward to the current GRASS 6 release mean that the interface is being rewritten, and this provides the opportunity to adapt it to the sp classes. Because GRASS provides the same kinds of data as sp classes handle, and relies on much of the same open source infrastructure (PROJ.4, GDAL, OGR), this step seems sensible. Work is still at an early stage, and some of the strucures used for GRASS 5 will not necessarily be retained. It does seem sensible, however, to try to respect the current region in GRASS to avoid handling raster data with different resolutions or extents. Like the GRASS 5 interface, R is assumed to be running within GRASS:
The metadata are accessed and available, but are not (yet) used to structure the sp class objects, here a SpatialGridDataFrame object filled with data from two spearfish layers. Here is a plot of the elevation data:
In addition, we can show what is going on inside the objects read into R:
List of 26
$ GISDBASE : chr "/home/rsb/topics/grassdata"
$ LOCATION_NAME: chr "spearfish57"
$ MAPSET : chr "rsb"
$ DEBUG : chr "0"
$ GRASS_GUI : chr "text"
$ projection : chr "1 (UTM)"
$ zone : chr "13"
$ datum : chr "nad27"
$ ellipsoid : chr "clark66"
$ north : num 4928010
$ south : num 4913700
$ west : num 589980
$ east : num 609000
$ top : num 1
$ bottom : num 0
$ nsres : num 30
$ nsres3 : num 30
$ ewres : num 30
$ ewres3 : num 30
$ tbres : num 1
$ rows : int 477
$ rows3 : int 477
$ cols : int 634
$ cols3 : int 634
$ depths : int 1
$ proj4 : chr "+proj=utm +zone=13 +a=6378206.4 +rf=294.9786982 +no_defs +nadgrids=/home/rsb/topics/grass61/grass-6.1.cvs/etc/nad/conus"
Object of class SpatialGridDataFrame
Coordinates:
min max
coords.x1 589980 609000
coords.x2 4913700 4928010
Is projected: TRUE
proj4string : [+proj=utm +zone=13 +a=6378206.4 +rf=294.9786982 +no_defs +nadgrids=/home/rsb/topics/grass61/grass-6.1.cvs/etc/nad/conus]
Number of points: 2
Grid attributes:
cellcentre.offset cellsize cells.dim
1 589995 30 634
2 4913715 30 477
Data attributes:
geology elevation.dem
sandstone:74959 Min. : 1066
limestone:61355 1st Qu.: 1200
shale :46423 Median : 1316
sand :36561 Mean : 1354
igneous :36534 3rd Qu.: 1488
(Other) :37636 Max. : 1840
NA's : 8950 NA's :10101
|
1 metamorphic 11693
2 transition 142
3 igneous 36534
4 sandstone 74959
5 limestone 61355
6 shale 46423
7 sandy shale 11266
8 claysand 14535
9 sand 36561
* no data 8950
Finally, a SpatialGridDataFrame object is written back to GRASS:
+----------------------------------------------------------------------------+
| Layer: sqdemSP Date: Sun May 14 21:59:26 2006 |
| Mapset: rsb Login of Creator: rsb |
| Location: spearfish57 |
| DataBase: /home/rsb/topics/grassdata |
| Title: ( sqdemSP ) |
|----------------------------------------------------------------------------|
| |
| Type of Map: raster Number of Categories: 255 |
| Data Type: FCELL |
| Rows: 477 |
| Columns: 634 |
| Total Cells: 302418 |
| Projection: UTM (zone 13) |
| N: 4928010 S: 4913700 Res: 30 |
| E: 609000 W: 589980 Res: 30 |
| Range of data: min = 32.649654 max = 42.895222 |
| |
| Data Source: |
| |
| |
| |
| Data Description: |
| generated by r.in.gdal |
| |
| |
+----------------------------------------------------------------------------+