CODEWITHSUNDEEP
rshapesspatialraster
Jian Zhang
Jian Zhang

Reputation: 1217

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)

Upvotes: 2

Views: 3627

Answers (1)

Jeffrey Evans
Jeffrey Evans

Reputation: 2417

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)

Upvotes: 6

Related Questions