When trying to rasterize a SpatialLine .shp, I get a raster with points instead of lines?

Question

EDIT: I've solved the problem, the resolution was simply not high enough to properly show the lines in the plot. My code now looks like this and works fine:

r <- raster(ncol=1000, nrow=1000, extent(sa))
res(r) <- 20
rivers <- rasterize(lines, r)

---

I've been trying (for two days now), to rasterize a SpatialLines .shp file in R, but every time I load the new raster, I get points instead of lines. I've tried gdal_rasterize and different values and fields, but nothing seems to work. This is a croped version of what it should look like with my DEM underneath:

But instead, it looks like this:

This is my code:

sa <- raster("studyarea.tif")
lines <- readOGR("Gewässer_Uckermark/gewnet25_bb_a.shp") 
r <- raster(ncol=1000, nrow=1000, extent(sa))
rivers <- rasterize(lines, r)
rivers[!(is.na(rivers))] <- 1

You can get the .shp file from here: https://geobroker.geobasis-bb.de/gbss.php?MODE=GetProductInformation&PRODUCTID=B9D461F1-99A1-4C10-97B4-9C36C0BD40B9

Answer 1

You could maybe use alternate rasterization package like velox (there is also fasterize but it is limited to polygons).

Be aware also that R raster plot displays a subset of pixels (typically a regular grid sample of 500000 pixels), so it may miss your tiny pixels (parameter maxpixels=500000 of plot for rasters).

Answer 2

This is caused by the plotting function which only plots a sample of the pixels for speed reasons. You can tweak it.

Here's a reproducible example taken from one feature from your shape: These coordinates are part of one feature in your data:

xy = structure(c(417476.8675, 417468.8832, 417473.037, 417458.7495, 
417352.3868, 417328.5742, 417318.52, 417303.7033, 417293.8735, 
5887737.1296, 5887834.9375, 5887938.579, 5888021.6583, 5888042.825, 
5888022.7166, 5887966.6249, 5887940.1665, 5887907.6158), .Dim = c(9L, 
2L), .Dimnames = list(NULL, c("x", "y")))

and we make a SpatialLines object out of it:

ltest = SpatialLines(list(Lines(list(Line(xy)),ID=1)))

we'll use a raster with a much larger extent:

e = extent(c(xmin=414500,xmax=424500,ymin=5882000,ymax=5890000))

Now try and rasterize over a fine grid on that raster:

r5 <- raster(ncol=4000, nrow=4000, extent(e))
rtest = rasterize(ltest,r5)
plot(ltest)
plot(rtest,add=TRUE)

That only shows a few pixels. But...

> sum(rtest[],na.rm=TRUE)
[1] 299

tells me there's 299 pixels. Change the maxpixels value to plot and:

> plot(ltest)
> plot(rtest,add=TRUE,maxpixels=ncell(rtest))

it looks fine:

Note the pixels seem much smaller than in the first plot, because I think when it subsamples them it also rescales them in size so they are visible in a full plot - that makes them way too big in a zoomed plot.

Any analysis you do with this raster will work fine, its only plotting that's affected and then you can fix that with maxpixels. Its not a bug despite what my earlier version may have said!