CODEWITHSUNDEEP
rshapefiler-raster
snoops
snoops

Reputation: 37

Add SpatialPointsDataFrame and SpatialLinesDataFrame to raster plot

I need to convert a shape file into a raster and have absolutely no clue where to even start.

If anyone could help me out I'd really appreciate!

Update: I have found out about the 'readOGR'-function, but everytime I use it I get the following message:

Warning messages:
1: In readOGR(dsn = "C:/Users/Giaco/Dropbox/Random Walk/Waterbodies",  :
  Dropping null geometries: 308, 309
2: In readOGR(dsn = "C:/Users/Giaco/Dropbox/Random Walk/Waterbodies",  :
  Z-dimension discarded

Can somebody tell me what that means?

Edit: 

altdata <- raster("altitude.tif")        
plot(altdata)
Lotic <- readOGR(dsn="C:/Users/Giaco/Dropbox/Random Walk/Waterbodies",layer="Lotic")
Lentic <- readOGR(dsn="C:/Users/Giaco/Dropbox/Random Walk/Waterbodies",layer="Lentic")

How can I plot the raster "altdata", the SpatialPointsDataFrame "Lentic" and the SpatialLinesDataFrame "Lotic" all in one plot ?

Edit: 

altdata
class       : RasterLayer 
dimensions  : 1286, 963, 1238418  (nrow, ncol, ncell)
resolution  : 15, 15  (x, y)
extent      : 90938.29, 105383.3, 190000, 209290  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=39.66666666666666 +lon_0=1 +k=1 +x_0=200000 +y_0=300000 +ellps=intl +pm=-9.131906111111112 +units=m +no_defs 
data source : C:\Users\Giaco\Dropbox\Random Walk\altitude.tif 
names       : altitude 
values      : -32768, 32767  (min, max)

> Lentic 
class       : SpatialPointsDataFrame 
features    : 182 
extent      : -108473.2, -95455.86, -107870.9, -91344.22  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=39.66825833333333 +lon_0=-8.133108333333334 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
variables   : 3
names       : Presence,     Type, Accessible 
min values  :        0, Fountain,          0 
max values  :        1,     Well,          2 
> Lotic
class       : SpatialLinesDataFrame 
features    : 317 
extent      : -108956.5, -93832.44, -108979.5, -90747.34  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=39.66825833333333 +lon_0=-8.133108333333334 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
variables   : 1
names       : Presence 
min values  :        0 
max values  :        1

Upvotes: 0

Views: 1418

Answers (1)

Jeffrey Evans
Jeffrey Evans

Reputation: 2397

Using base plot, there are two way to overlay vector data on a raster. First plot the raster, then you can either call plot for each feature class using the add=TRUE argument. Alternately, you can use the points and lines functions which will also add to the current plot.

Create some example data

library(raster)
library(sp)

 r <- raster(nrows=180, ncols=360, xmn=571823.6, xmx=616763.6, ymn=4423540, 
             ymx=4453690, resolution=270, crs = CRS("+proj=utm +zone=12 +datum=NAD83 
             +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
 r[] <- runif(ncell(r))
 pts <- sampleRandom(r, 10, na.rm = TRUE, sp = TRUE)

Plot the raster then call plot again, with add=TRUE, to add points. 

plot(r)
  plot(pts, pch=20, cex=2, col="red", add=TRUE)

Or, plot the raster and use points to add the point feature class. 

plot(r)
  points(pts, pch=20, cex=2, col="red")

Edit: Your extents between the raster and vector feature classes do not overlap.

We can create SpatialPolygons using the extent of your objects and an example raster (with a uniform value of 1). 

library(raster)
proj <- sp::CRS("+proj=tmerc +lat_0=39.66666666666666 +lon_0=1 +k=1 +x_0=200000 +y_0=300000  
                 +ellps=intl +pm=-9.131906111111112 +units=m +no_defs")
r.ext <- as(extent(90938.29, 105383.3, 190000, 209290), "SpatialPolygons")
  proj4string(r.ext) <- proj
pt.ext <- as(extent(-108473.2, -95455.86, -107870.9, -91344.22), "SpatialPolygons")
  proj4string(pt.ext) <- proj
line.ext <- as(extent(-108956.5, -93832.44, -108979.5, -90747.34),  "SpatialPolygons")
    proj4string(line.ext) <- proj
r <- raster(r.ext, resolution = c(15,15), crs = proj)
  r[] <- rep(1, ncell(r))

Here we see that if we plot the raster and then the point and line extent polygons, you cannot see them.

plot(r)
  plot(pt.ext, add=TRUE)
  plot(line.ext, add=TRUE)

However, if we plot the line and point extent polygons they overlay just fine.

plot(line.ext)
  plot(pt.ext,add=TRUE)

If we limit the raster extent to the extent of the line object, you should see the raster sub-region, but do not. And, if you try to crop the raster you will receive an "Error in .local(x, y, ...) : extents do not overlap" error. 

plot(line.ext)
  plot(r, add=TRUE)

Upvotes: 2

Related Questions