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:

enter image description here

But instead, it looks like this: enter image description here

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


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).

| improve this answer | |
  • Aaagh. I knew I'd seen this before. I've augmented my answer to add this. – Spacedman May 11 at 10:42

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)

enter image description here

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:

enter image description here

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!

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.