Software
From spatial-analyst.net
his article provides some basic instructions on how to install and do first steps in software the R + Open Source Desktop GIS (osgeo) + Google Earth bundle (this article was specifically written for the MS Windows OS users) that can be used to run analysis of various spatio-temporal data. This bundle allows full GIS+statistics integration and can support practically about 80% of processing/visualization capabilities available in commercial packages such as ArcInfo/Map or Idrisi (compare the web-traffic trends for ESRI and R+osgeo-project). The advantage of combining R with open source GIS is that you will be able to process and visualize even large datasets (R is not really fit for processing of large data). In addition, packages such as ILWIS and SAGA allow you to input, edit and visually explore geographical data, before and after the actual statistical analysis. Note that there is a large number of open GIS alternatives and commercial GIS packages that could also be linked to R with equal success. For example, GRASS GIS (versions >6.3 now also available for MS Windows machines) is another good candidate for GIS+statistics integration. Within GRASS, you can load R and then use its functionality on the GRASS map layers, or run R on top of GRASS (with the help of the spgrass6 package). Another alternative for GIS+R integration is the QGIS, which has among its main characteristics a python console, and a very elaborated way of adding Python plugins, which is already used for an efficient R plugin (manageR).
Contents |
R
Installation
To install R under Windows, download and run an installation exe from the R project homepage. This will install R for window with a GUI. After you start R, you will first need to set-up the working directory and install additional packages. To run geostatistical analysis in R, you will need to add the following R packages: gstat (gtat in R), rgdal (GDAL import of GIS layers in R), sp (operations on maps in R), foreign, spatstat (spatial statistics in R) and maptools. To install these packages you should do the following. First start the R GUI, then select Packages > Load package from the main menu. Note that, if you wish to install a package on the fly, you will need to select a suitable CRAN mirror from where it will download and unpack a package. You can also install a number of packages used for the spatial analysis by typing e.g.:
> install.packages(c("rgdal","maptools","gstat","adehabitat","akima","spatstat","RSAGA"))
Another quick way to get all packages used in R to do spatial analysis (as explained in the ASDAR book by Bivand et al. 2008) is to install the ctv package and then execute the command:
> install.packages("ctv")
> library(ctv)
> install.views("Spatial")
Although a package is installed, you will still need to load it into your workspace (every time you start R) before you can use its functionality. A package is commonly loaded using e.g.:
> library(gstat)
First steps in R
R is today identified as one of the fastest growing and most comprehensive statistical computing tools/communities. It practically offers statistical analysis and visualisation of unlimited sophistication. A user is not restricted to a small set of procedures or options, and because of the contributed packages, users are not limited to one method of accomplishing a given computation or graphical presentation (Bivand et al., 2008; Rossiter, 2007; Murrell, 2006). R became attractive for analysis of spatial data mainly due to the recent integration of the geostatistical tools (e.g. gstat, geoR) and tools that allow R computations with spatial data layers (e.g. sp/maptools, rgdal, r-gis). A complete listing of all R packages that can be used to run analysis on spatial data in R is maintained by Roger Bivand. The 'Spatial' packages can be nicely combined with e.g. the Environmentrics packages. The interactive graphics in R is also increasingly powerful.
Although there are few GUI's being developed to run analysis in R (e.g. Rcmdr), R is (and will stay) a command-based software. So, the first step you should consider is to get some initial training on R syntax. You may consult some literature first, but also try to follow 1-2 days training courses. The second useful step to get into R is to try to obtain a more user-friendly R scripting edition such as as TINN-R or JGR. You can use this Rprofile.site (copy to your "R/etc/" directory) and this Rconsole setting to produce the same GUI configuration as shown below.
If a single argument in the command is incorrect, inappropriate or mistyped, you will get an error message. If the error message is not helpful, then try receiving more help about some operation by typing e.g.:
> help(krige)
or
> ?krige
commands or via the Html help files. You may also wish to register to the special interest groups such as R-sig-Geo or similar and subscribe to their mailing lists. Gstat-info is the mailing list of the gstat package and also offers many interesting information.
These are some easy to follow R tutorials:
- Simple R
- Quick R
- Introduction to R by D.G. Rossiter
- Introductory course on Spatial Data Analysis with R by R. Bivand
A brief guide to learn how to make scripts in R is available in this article. But please, consider also purchasing some of the official R books and/or attending some of the short courses on R.
To download the file (DEM25m) used in this exercise, we can do:
> download.file("http://spatial-analyst.net/DATA/DEM25m.zip", destfile=paste(getwd(),"DEM25m.zip",sep="/"))
> fname <- zip.file.extract(file="DEM25m.asc", zipname="DEM25m.zip")
> file.copy(fname, "./DEM25m.asc", overwrite=TRUE)
Next, we can import the ArcInfo ASCII map to R using the rgdal package:
> library(rgdal)
> dem25m <- readGDAL("DEM25m.asc")
After we attach the correct coordinate system string to this map, we can export it to e.g. ILWIS GIS by using:
> proj4string(dem25m) <- CRS("+proj=tmerc +lat_0=0 +lon_0=18 +k=0.9999 +x_0=6500000 +y_0=0 +ellps=bessel
+ +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824 +units=m")
> writeGDAL(dem25m[1], "dem25m.mpr", "ILWIS")
A serious limitation of R, considering GIS data, is the memory handling issue. You will soon discover that it is almost impossible to visualize and process raster maps which are >1--5 million pixels. In addition, R has still limited capabilities for visual data exploration, and almost no capabilities for data input/editing. You can overcome these problems by combining R with the open source GIS, such as SAGA, GRASS or ILWIS.
ILWIS GIS
Installing ILWIS
ILWIS is an acronym for Integrated Land and Water Information System, a stand alone integrated GIS package developed at the International Institute of Geoinformation Science and Earth Observations (ITC), Enschede, Netherlands. ILWIS was originally built for educational purposes and low-cost applications in developing countries. Its development started in 1984 and the first version (DOS version 1.0) was released in 1988. ILWIS 2.0 for Windows was released at the end of 1996, and a more compact and stable version 3.0 (WIN 95) was released by mid 2001. From 2004, ILWIS was distributed solely by ITC as shareware at a nominal price. From July 2007, ILWIS shifted to open source and ITC will not provide support for its further development. Further development of ILWIS can be mainly followed via the 52 North initiative and via the ILWIS mailing list.
ILWIS (3.5) offers a range of image processing, vector, raster, geostatistical, statistical, database and similar operations. In addition, a user can create new scripts, adjust the operation menus and even build Visual Basic, Delphi, or C++ applications that will run at top of ILWIS and use its internal functions. The real advantage of ILWIS is that all operations can be run from command line using the ILWIS syntax; list of commands can then be combined in ILWIS scripts to automate data processing. To learn more about ILWIS, its design, operations and syntax, try to obtain some of the following literature:
- Hengl, T., Maathuis, B.H.P., Wang, L., 2008. Geomorphometry in ILWIS. In: Hengl, T. and Reuter, H.I. (Eds), Geomorphometry: Geomorphometry: Concepts, Software, Applications. Developments in Soil Science, vol. 33, Elsevier, 309-330 pp.
- Unit Geo Software Development, 2001. ILWIS 3.0 Academic User’s Guide. International Institute for Geo-Information Science and Earth Observation (ITC), Enschede, 530 pp.
- ILWIS Help documentation.
The most recent version of ILWIS -- 3.5 (an installation can be obtained from here) -- is developed as a MS Visual 2008 project. ILWIS user interface and ILWIS analytical functionality have now been completely separated making it easier to write server side applications for ILWIS. This allows us to control ILWIS also from R, e.g. by setting the location of ILWIS on your machine:
ILWIS <- "C:\\Progra~1\\N52\\Ilwis35\\IlwisClient.exe -C"
ILWIS syntax
ILWIS has its original syntax to run spatial analysis. In many cases, it is fairly simple to understand and reproduce the syntax. E.g. to resample a raster map from one to another grid system, we can can use the MapResample command:
dem25m_ll_c.mpr = MapResample(dem25m_c.mpr, geoarc, BiLinear)
This will create a new map using a Bilinear resampling and geoarc as the new grid definition. To combine ILWIS and R commands in R, we use:
> shell(cmd=paste(ILWIS, " crgrf geoarc.grf ",geoarc@grid@cells.dim[[2]],"
+ ",geoarc@grid@cells.dim[[1]]," -crdsys=LatlonWGS84 -lowleft=(",geoarc@grid@cellcentre.offset[[1]],"
+ ,",geoarc@grid@cellcentre.offset[[2]],") -pixsize=",geoarc@grid@cellsize[[1]],sep=""),
+ wait=F)
> shell(cmd=paste(ILWIS, "dem25m_ll_c.mpr = MapResample(dem25m_c.mpr, geoarc, BiLinear)"), wait=F)
> shell(cmd=paste(ILWIS, "open dem25m_ll_c.mpr -noask"), wait=F)
As a result, you will get a screen as shown on the right.
The advantage of using ILWIS (as compared to the spplot command of the sp package in R) is that it allows direct visual exploration and editing of spatial data. It can also handle relatively large grids, which might be more complicated to visualize in R. Note also that ILWIS extends the image processing functionality (especially image filtering, resampling, geomorphometric analysis and similar) of R that is, at the moment, limited to only few experimental packages (e.g. biOps, rimage; raster, see also the list of r-gis projects under development).
SAGA GIS
Installing SAGA
SAGA stands for System for Automated Geoscientific Analyses is a full-fledged GIS, and many of its features have some relation with geomorphometry, which makes it an ideal tool for operational work, but also for GIS training purposes. It is an open-source GIS with support for raster and vector data. It includes a large set of geoscientific algorithms, and is especially powerful for the analysis of DEMs.
SAGA has been under development since 2001 at the University of Göttingen (the SAGA development team, including the key developer Olaf Conrad, has since moved to University of Hamburg), Germany, with aim of simplifying the implementation of new algorithms for spatial data analysis within a framework that immediately allows their operational application. Therefore, SAGA targets not only the pure user but also the developer of geo-scientific methods. SAGA has its roots in DiGeM, a small program specially designed for the extraction of hydrological land-surface parameters. In 2004 most of SAGA’s source code was published using an Open Source Software (OSS) license. To learn more about SAGA, try obtaining some of the following references:
- Böhner, J., McCloy, K.R., Strobl, J. (Eds.), 2006. SAGA — Analysis and Modelling Applications. Göttinger Geographische Abhandlungen, Heft 115. Verlag Erich Goltze GmbH, Göttingen, 117 pp.
- Brenning, A. 2008. Statistical geocomputing combining R and SAGA: The example of landslide susceptibility analysis with generalized additive models. In: J. Böhner, T. Blaschke & L. Montanarella (eds.), SAGA - Seconds Out (Hamburger Beiträge zur Physischen Geographie und Landschaftsökologie, 19), 23-32.
- Conrad, O. 2007. SAGA - Entwurf, Funktionsumfang und Anwendung eines Systems für Automatisierte Geowissenschaftliche Analysen. electronic doctoral dissertation, University of Göttingen.
- Olaya, V., Conrad, O., 2008. Geomorphometry in SAGA. In: Hengl, T. and Reuter, H.I. (Eds), Geomorphometry: Concepts, Software, Applications. Developments in Soil Science, vol. 33, Elsevier, 293-308 pp.
SAGA can be obtained from a Source Forge repository; the most recent version is available both for windows and linux operating systems. A list of available libraries and modules can be also browsed from the SAGA GIS homepage.
First steps in SAGA
First, you need to load the RSAGA library by:
> library(RSAGA)
Note that loading a library you will not be able to run any of the SAGA operations if also SAGA binaries and libraries are not installed on the same machine (both SAGA and ILWIS are external GIS packages and NOT built-in R packages). You can at any time set the location of your SAGA environment by typing:
> rsaga.env(path="C:/Progra~1/saga_vc")
To find out what a certain module (the best way to browse all possible modules is to open SAGA and look at the left menu tab that shows all module libraries available) in SAGA does, you can execute:
> rsaga.get.modules("grid_gridding")
which will show this:
$grid_gridding code name interactive 1 0 Inverse Distance FALSE 2 1 Nearest Neighbour FALSE 3 2 Modifed Quadratic Shepard FALSE 4 3 Shapes to Grid FALSE 5 4 Triangulation FALSE 6 NA <NA> FALSE 7 NA <NA> FALSE Warning message: NAs introduced by coercion
To get further usage of the 3rd function of the grid_gridding module, you can execute:
> rsaga.get.usage("grid_gridding", 2),
which then gives the complete list of parameters (obligatory and optional) that need to be set before this command can be executed:
SAGA CMD 2.0.3
library path: C:/Progra~1/saga_vc/modules
library name: grid_gridding
module name : Shapes to Grid
Usage: 3 [-GRID <str>] -INPUT <str> [-FIELD <num>] ...
-GRID:<str> Grid
Data Object (optional output)
-INPUT:<str> Shapes
Shapes (input)
...
Once you have learned about a specific function, you can run it via the rsaga.geoprocessor. For example to extract a stream network from a DEM, we can run (note that we first convert the ArcInfo ASCII DEM to the SAGA format):
> rsaga.esri.to.sgrd(in.grids="dem25m.asc", out.sgrds="dem25m.sgrd", in.path=getwd())
> rsaga.geoprocessor(lib="ta_channels", module=0, param=list(ELEVATION="dem25m_f.sgrd",
+ CHNLNTWRK="channel_ntwrk.sgrd", CHNLROUTE="channel_route.sgrd", SHAPES="channels.shp",
+ INIT_GRID="dem25m_f.sgrd", DIV_CELLS=3, MINLEN=40))
> channels <- readOGR("channels.shp", "channels")
> spplot(channels["LENGTH"], col.regions=bpy.colors())
which will produce a screen as shown on the right (above). Note that the rsaga.geoprocessor only serves to pass the parameters to SAGA command line (saga_cmd.exe). If you are loading a list of grids to some SAGA process, they need to be written as a list:
# export grids back to ArcInfo ASCII format:
> rsaga.sgrd.to.esri(in.sgrds=c("channel_ntwrk.sgrd","channel_route.sgrd"),
+ out.grids=c("channel_ntwrk.asc","channel_route.asc"), out.path=getwd(), prec=1)
the rsaga.geoprocessor would expect grids separated by ";", e.g. GRIDS="DSM1.sgrd;DSM2.sgrd;DSM3.sgrd". Sometimes you might receive an error message or computation will fail, which possibly means a bug in SAGA GIS. Such problems you should report directly via the SAGA user's forum (and not via the R-sig-geo mailing list). Read more on how to extract land-surface parameters and automate analysis of DEMs using SAGA.
GRASS GIS
Installing GRASS (MS Windows)
GRASS GIS, now one of the eight initial Software Projects of the Open Source Geospatial Foundation, is probably the most known open source GIS software and does need to be specially introduced. Its functionality and usage are described in detail in Neteler et al. (2008); a wiki-based tutorial is available at http://grass.osgeo.org/wiki/. GRASS itself is a collection of environment variables (they vary from version to version), which are populated when GRASS is started --- that is a GRASS shell is started within the running command shell. GRASS is a much larger projected than SAGA considering the number of developers/institutions involved, although their functionality considering the DEM analysis is about similar. The latest version of GRASS (6.3.) contains over 350 routines.
Although originally a Linux-based project, the most recent version of GRASS (6.3; development version) is now also available for MS Windows machines. The coming version 6.4. is expected to be even more stable under Windows. To install GRASS, simply obtain the installation exe and install GRASS to some default location e.g. "C:/GRASS/"
Controlling GRASS from R
The GRASS+R integration requires the spgrass6 package. This needs to be installed with all necessary dependencies:
> install.packages("spgrass6", dependencies=T)
> library(spgrass6)
Next, you need to set the location of your GRASS installation:
> loc <- initGRASS("C:/GRASS", home=tempdir())
> loc
gisdbase c:/Temp
location file678418be
mapset file3d6c4ae1
rows 1
columns 1
north 1
south 0
west 0
east 1
nsres 1
ewres 1
projection NA
This will automatically generate a temporary GRASS gisdbase (and create a temporary file). Note that this database contains no maps. These can be imported using:
> execGRASS("r.in.gdal", flags="o", parameters=list(input="DEM25m.asc", output="DEM"))
WARNING: Over-riding projection check
RINGDA~1 complete.
> execGRASS("g.region", parameters=list(rast="DEM"))
> gmeta6()
gisdbase c:/Temp
location file678418be
mapset file3d6c4ae1
rows 149
columns 147
north 5074275
south 5070550
west 6551872
east 6555547
nsres 25
ewres 25
projection NA
To get a description of some GRASS function and parameters connected, you can use e.g.:
> parseGRASS("r.in.gdal")
Once you have imported a map to GRASS format, you can also open and visualize this DEM using the GRASS GUI. In the next step, we can extract the drainage network:
> execGRASS("r.watershed", flags="overwrite",
+ parameters=list(elevation="DEM", stream="stream", threshold=as.integer(50)))
SECTION 1a (of 5): Initiating Memory.
SECTION 1b (of 5): Determining Offmap Flow.
SECTION 2: A * Search.
SECTION 3: Accumulating Surface Flow.
SECTION 4: Watershed determination.
SECTION 5: Closing Maps.
# thin the raster map so it can be converted to vectors:
> execGRASS("r.thin", parameters=list(input="stream", output="streamt"))
File stream -- 149 rows X 147 columns
Bounding box: l = 2, r = 148, t = 2, b = 150
Pass number 1
Deleted 105 pixels
Pass number 2
Deleted 0 pixels
Thinning completed successfully.
Output file 149 rows X 147 columns
Window 149 rows X 147 columns
# convert to vectors:
> execGRASS("r.to.vect", parameters=list(input="streamt", output="streamt", feature="line"))
WARNING: Default driver / database set to:
driver: dbf
database: $GISDBASE/$LOCATION_NAME/$MAPSET/dbf/
Extracting lines...
Building topology ...
Registering lines:
287 primitives registered
Building areas:
0 areas built
0 isles built
Attaching islands:
Attaching centroids:
Topology was built.
Number of nodes : 266
Number of primitives: 287
Number of points : 0
Number of lines : 287
Number of boundaries: 0
Number of centroids : 0
Number of areas : 0
Number of isles : 0
RTOVEC~1 complete.
Finally, we can import the GRASS maps to R using the spgrass6 package:
> streamt <- readVECT6("streamt")
Exporting 287 points/lines...
287 features written
OGR data source with driver: ESRI Shapefile
Source: "c:/Temp/file678418be/file3d6c4ae1/.tmp", layer: "streamt"
with 287 rows and 3 columns
Feature type: wkbLineString with 2 dimensions
> plot(streamt)
|
