Reading and writing spatial data¶
Introduction¶
Reading and writing spatial data is complicated by the fact that there are many different file formats. However, there are a few formats that are most common that we discuss here.
Vector files¶
The shapefile is the most commonly used file format for vector data.
It is trivial to read and write such files. Here we use a shapefile that
comes with the raster package.
Reading¶
We use the system.file function to get the full path name of the
file’s location. We need to do this as the location of this file depends
on where the raster package is installed. You should not use the
system.file function for your own files. It only serves for creating
examples with data that ships with R.
library(raster)
filename <- system.file("external/lux.shp", package="raster")
filename
## [1] "C:/soft/R/R-devel/library/raster/external/lux.shp"
Now we have the filename we need we use the shapefile function. This
function comes with the raster package. For it to work you must also
have the rgdal package.
s <- shapefile(filename)
s
## class : SpatialPolygonsDataFrame
## features : 12
## extent : 5.74414, 6.528252, 49.44781, 50.18162 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 5
## names : ID_1, NAME_1, ID_2, NAME_2, AREA
## min values : 1, Diekirch, 1, Capellen, 76
## max values : 3, Luxembourg, 12, Wiltz, 312
The shapefile function returns Spatial*DataFrame objects. In
this case a SpatialPolygonsDataFrame. It is important to recognise
the difference between this type of R object
(SpatialPolygonsDataFrame), and the file (shapefile) that was used to
create it.
For other formats, you can use readOGR function in package
rgdal.
Writing¶
You can also write shapefiles using the shapefile method. In stead
of a filename, you need to provide a vector type Spatial* object as
first argument and a new filename as a second argument. You can add
argument overwrite=TRUE if you want to overwrite an existing file.
outfile <- 'test.shp'
shapefile(s, outfile, overwrite=TRUE)
For other formats, you can use writeOGR function in package
rgdal.
Raster files¶
The raster package can read and write several raster file formats.
Reading¶
Again we need to get a filename for an example file.
f <- system.file("external/rlogo.grd", package="raster")
f
## [1] "C:/soft/R/R-devel/library/raster/external/rlogo.grd"
Now we can do
r1 <- raster(f)
r1
## class : RasterLayer
## band : 1 (of 3 bands)
## dimensions : 77, 101, 7777 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 0, 101, 0, 77 (xmin, xmax, ymin, ymax)
## crs : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## source : rlogo.grd
## names : red
## values : 0, 255 (min, max)
Note that r1 is a RasterLayer of the first “band” (layer) in the
file (out of three bands (layers)). We can request another layer.
r2 <- raster(f, band=2)
r2
## class : RasterLayer
## band : 2 (of 3 bands)
## dimensions : 77, 101, 7777 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 0, 101, 0, 77 (xmin, xmax, ymin, ymax)
## crs : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## source : rlogo.grd
## names : green
## values : 0, 255 (min, max)
More commonly, you would want all layers in a single object. For that
you can use the brick function.
b <- brick(f)
b
## class : RasterBrick
## dimensions : 77, 101, 7777, 3 (nrow, ncol, ncell, nlayers)
## resolution : 1, 1 (x, y)
## extent : 0, 101, 0, 77 (xmin, xmax, ymin, ymax)
## crs : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## source : rlogo.grd
## names : red, green, blue
## min values : 0, 0, 0
## max values : 255, 255, 255
Or you can use stack, but that is less efficient in most cases.
s <- stack(f)
s
## class : RasterStack
## dimensions : 77, 101, 7777, 3 (nrow, ncol, ncell, nlayers)
## resolution : 1, 1 (x, y)
## extent : 0, 101, 0, 77 (xmin, xmax, ymin, ymax)
## crs : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## names : red, green, blue
## min values : 0, 0, 0
## max values : 255, 255, 255
The same approach holds for other raster file formats, including GeoTiff, NetCDF, Imagine, and ESRI Grid formats.
Writing¶
Use writeRaster to write raster data. You must provide a Raster*
object and a filename. The file format will be guessed from the filename
extension (if that does not work you can provide an argument like
format=GTIFF). Note the argument overwrite=TRUE and see
?writeRaster for more arguments, such as datatype= to set the
datatype (e.g., integer, float).
x <- writeRaster(s, 'output.tif', overwrite=TRUE)
x
## class : RasterBrick
## dimensions : 77, 101, 7777, 3 (nrow, ncol, ncell, nlayers)
## resolution : 1, 1 (x, y)
## extent : 0, 101, 0, 77 (xmin, xmax, ymin, ymax)
## crs : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## source : output.tif
## names : output.1, output.2, output.3
## min values : 0, 0, 0
## max values : 255, 255, 255