Search code examples
rshapesspatialraster

Convert GIS shape file to raster in R?


I need to convert shapefiles into raster format.

I used the function "rasterize" in R package "raster", but the result does not look correct.

tst <- rasterize(shpfile, r, fun="count")  
Found 5 region(s) and 5 polygon(s)

There is no gird with occurrence records:

range(tst[],na.rm=TRUE)
[1]  Inf -Inf
Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf

sum(tst[],na.rm=TRUE)
[1] 0

The R script that I wrote:

# download the GIS shape file
download.file("http://esp.cr.usgs.gov/data/little/abiebrac.zip",
              destfile = "abiebrac.zip")

# unzip
unzip("abiebrac.zip")

# Read shapefile into R
library(maptools)

shpfile <- readShapePoly("abiebrac")
extent(shpfile)

# mapping
plot(shpfile)

library(maps)
map("world",xlim=c(-180,-50),ylim=c(7,83),add=FALSE)
plot(shpfile,add=TRUE,lwd=10)


# rasterize
library(raster)

GridSize <- 0.5 # degree
r <- raster(ncol= round(abs(-180-(-50))/GridSize), 
            nrow=round(abs(83-7)/GridSize))
extent(r) <- extent(c(-180, -50, 7, 83))
tst <- rasterize(shpfile, r, fun="count")

# summary
sum(tst[],na.rm=TRUE)
range(tst[],na.rm=TRUE)

# mapping
plot(tst,col="red",main="abiebrac")  
map("world",xlim=c(-180,-50),ylim=c(7,83),add=TRUE)

Solution

  • I am not sure why you are using "count" in the fun argument but in this case, because there is no overlap, it is producing NA results. You also need to define an attribute field in the spatialPolygonDataFrame object to assign values to your raster. You can also pull the extent directly from the sp object.

    This code seems to yield what you want.

    require(raster)
    require(rgdal)
    require(sp)
    
    setwd("D:/TMP")
    shpfile <- readOGR(getwd(), "abiebrac")
    r <- raster(extent(shpfile))
      res(r)=0.05
        r <- rasterize(shpfile, field="ABIEBRAC_", r)
    
    plot(r)
      plot(shpfile,lwd=10,add=TRUE)