I have the following raster (shp file using readORG):

structure(list(class = structure("SpatialPolygonsDataFrame", package = "sp"), 
    bbox = structure(c(4346000, 3819000, 4445000, 3867000), .Dim = c(2L, 
    2L), .Dimnames = list(c("x", "y"), c("min", "max"))), is.projected = TRUE, 
    proj4string = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs", 
    data = structure(c("Min.   :0  ", "1st Qu.:0  ", "Median :0  ", 
    "Mean   :0  ", "3rd Qu.:0  ", "Max.   :0  "), .Dim = c(6L, 
    1L), .Dimnames = list(c("", "", "", "", "", ""), "      Id"), class = "table")), .Names = c("class", 
"bbox", "is.projected", "proj4string", "data"), class = "summary.Spatial")

as well as the shp file containing lines:

structure(list(class = structure("SpatialLinesDataFrame", package = "sp"), 
    bbox = structure(c(4329488.96922647, 3429159.10800761, 4998503.48859431, 
    4055688.10547651), .Dim = c(2L, 2L), .Dimnames = list(c("x", 
    "y"), c("min", "max"))), is.projected = TRUE, proj4string = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs", 
    data = structure(c("Min.   :0  ", "1st Qu.:0  ", "Median :0  ", 
    "Mean   :0  ", "3rd Qu.:0  ", "Max.   :0  ", "Min.   :  448  ", 
    "1st Qu.:13229  ", "Median :28235  ", "Mean   :27205  ", 
    "3rd Qu.:40724  ", "Max.   :50608  ", "Min.   :  485  ", 
    "1st Qu.:13731  ", "Median :29399  ", "Mean   :28635  ", 
    "3rd Qu.:43159  ", "Max.   :53607  ", "Min.   :6519754  ", 
    "1st Qu.:8741416  ", "Median :9167928  ", "Mean   :8894830  ", 
    "3rd Qu.:9414708  ", "Max.   :9762259  ", "Fishing:121  ", 
    NA, NA, NA, NA, NA), .Dim = c(6L, 5L), .Dimnames = list(c("", 
    "", "", "", "", ""), c("      Id", "   trip_id", "  new_tr_id", 
    "     species", "  HELCOM_Gro")), class = "table")), .Names = c("class", 
"bbox", "is.projected", "proj4string", "data"), class = "summary.Spatial")

In this shp file with the lines, each lines has a unique identification number called new_tr_id.

My goal is to create a raster file with the number if lines crossing each cell. A line (so the same new_tr_id) can cross the same cell X times, it would be counted X times. I added the figure below as example:

In red are the values for each cell. I added the new_tr_id for each line to understand the counting method. For example in the cell surrounded with the blue color, the value is 2 because the lines with the new_tr_id 3 is crossing 2 times

I have some interesting discussions using the package raster but I still can not find a proper answer to this issue.

Some news: I tried the method described by @Spacedman below but it is cutting the lines crossing together: enter image description here

I don´t have the same result as @Spacedman but I am not that far anymore :)

  • "A line (so the same new_tr_id) can cross the same cell X times, it would be counted X times." Can you upload a figure showing that? – aldo_tapia Dec 1 '17 at 11:46
  • @aldo_tapia, just did it! – Floni Dec 1 '17 at 12:01
  • I think this task could be achieved converting raster to polygons using new_tr_id as value field... R interprets rasters as matrices, so could be difficult to do (and I don't know how to do it) – aldo_tapia Dec 1 '17 at 12:15
  • rasterize(lines, grid, fun = "count") will do it but you'll need a raster version of your grid, something like rasterFromXYZ(coordinates(summary_grid)) – mdsumner Dec 1 '17 at 12:19
  • Thanks @mdsumner. So ´rasterFromXYZ()´ is changing the shape file into raster version, right? – Floni Dec 1 '17 at 12:24

Here's a method using sf packages. So install sf and lwgeom and dplyr first.

polys is read from a polygon shapefile using polys = st_read("polys.shp") - lines is read from a line shapefile using lines=st_read("lines.shp")/

My lines has a unique FID field and my polygons have a unique id field.

First compute the intersection:

> pli = st_intersection(polys,lines)

Each row of this has the id and FID field of the relevant polygon and line. But in the cases where multiple line segments from the same line appear twice or more in the polygon return as MULTILINESTRING segments. For these we need to count the number of LINESTRING segments. I couldn't find a way to do this better than this:

countsegs =
    if("LINESTRING" %in% class(g)){

which returns 1 for a line string and the number of linestrings for a multilinestring. We continue:

> pli$segs = sapply(pli$geometry, countsegs)

Now there's a column with the number of segments in each polygon. Use some dplyr package pipes and whatnot to sum over polygon id to get a data frame of counts of segments in polygons. There's probably an easier way to strip the geometry out but whatevs:

> summed = pli %>% group_by(id) %>%
       summarise(total = sum(segs)) %>% select(id, total)
> summed = as.data.frame(summed)
> summed$geometry=NULL

Now join the totals onto the polygon spatial data frame:

> polys = left_join(polys, summed, c("id"="id"))
> plot(polys[,"total"])
> plot(lines$geom, add=TRUE)

Does this look right?

enter image description here

| improve this answer | |
  • The sf package, in conjuction with packages like dplyr, are really making this type of analyses quite elegant. Thank you for showing an example that is moving us towards this spatial object class. Do you think that it is time to start retooling our packages to use the sf class or is sp still prefered? – Jeffrey Evans Dec 1 '17 at 18:06
  • All you need is sf and as(thing,"Spatial") if you need sp classes for packages that work with them. – Spacedman Dec 1 '17 at 18:10
  • Thank you @spacedman. I will try ASAP. I can see that you actually cut each lines for each cells. I hope it can work for a large grid and dataset. We will know soon :) – Floni Dec 3 '17 at 10:42
  • thanks again @Spacedman, but it does not work with my sample. Is it possible to have your sample to try out and to do a match with mine? – Floni Dec 12 '17 at 7:37
  • @spacedman, it seems to work but the lines are also splitted when they cross each other. For example in the cells in yellow in your results, the value is not 4 but more! – Floni Dec 14 '17 at 7:24

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.