Generating a raster layer with intersecting line lengths as cell values

Question

I have a polyline vector layer and a template raster dataset. For this polyline vector layer I would like to generate a new raster dataset based on the template raster dataset, but in which each cell contains a value corresponding to the length of the line segment that crosses that cell. All other cells should contain a value Nodata. Heres a stupid example I made in Paint:

I am using python based tools including rasterio, geopandas, shapely and rasterstats. I would like to avoid using lower level libraries such as gdal and ogr if I can help it, but I will resort to them if it makes the difference between life and death.

What I am currently doing is the following:

• Extract template raster dataset by mask using polyline dataset
• reassign all values in extracted dataset a unique value

• polygonize raster dataset based on following function

  def polygonize_raster(dataset):
  
      # Read the dataset's valid data mask as a ndarray. Dataset is a rasterio read object open for reading
      mask = dataset.dataset_mask()
  
      array = dataset.read(1)
      generator = rasterio.features.shapes(source=array, mask=mask, transform=dataset.transform)
      # Extract feature shapes and values from the array
      geom_list = []
      for geom, value in generator:
          # Print GeoJSON shapes to stdout
          geom = shapely.geometry.shape(geom)
          geom_list.append(geom)
      return geom_list
  

• generate a new list with line lengths intersecting each polygon based on the following functions:

  def calc_length_in_cell(polygon_geom, linestring_geom):
  
      length = linestring_geom.intersection(polygon_geom).length
      return length
  
  
  def calc_length_in_cell_list(polygon_geom_list, multilinestring):
  
      #loops through all polygon geometries, calculates cell length inside and returns list with all lengths
      lengths = []
      for i,polygon_geom in enumerate(polygon_geom_list):
          lengths.append(calc_length_in_cell(polygon_geom, multilinestring))
      return lengths
  

• add the list of lengths to a geodataframe containing the polygonized raster (little boxes corresponding to each template raster cell).

• burn the above geodataframe into a new raster with the line lengths as burn value

I have 3 issues with this workflow:

• Smaller tests with square rasters containing about 9 cells with unique values seem to polygonize just fine, and each cell gets its own polygon. However for large datasets this stops being the case, as evinced by the image below despite the fact that adjacent cells have different values

• the calc_length_in_cell_list function above is a killer. It really blows up for large datasets since it has to perform an intersection for every polygon. Heavy stuff.

• the polygonized raster is technically invalid geometry as can be determined by this question.

I would like to either 1) come up with a different workflow or 2) tweak the current one such that this can be accomplished for large datasets with ease, and that the resulting raster with intersecting polyline lengths is correct.

Answer 1

One way to approximate this is by:

1. rasterizing to a resolution finer than the desired output
2. Aggregating the raster to the desired resolution

See the example below that first rasterizes the LineString to a 0.1 degrees resolution raster/array and then aggregates it with sum() to 1 degree

from rasterio import features
from affine import Affine
import numpy as np

# LineString going through France, Spain, Mediterranean sea
geom = {'coordinates': [[-0.703, 48.283],
                        [2.285, 47.576],
                        [0.966, 46.8],
                        [2.021, 46.195],
                        [1.318, 44.995],
                        [3.032, 44.308],
                        [4.658, 43.357],
                        [5.888, 40.946],
                        [-3.251, 40.913],
                        [11.865, 37.3],
                        [9.272, 39.943]],
        'type': 'LineString'}

# Define affine with 0.1 degrees resolution
aff = Affine(0.1, 0, -5,
             0, -0.1, 50)

# Rasterize at 0.1 degrees resolution
arr_0 = features.rasterize([geom], out_shape=(200, 200),
                           fill=0, transform=aff,
                           dtype=np.int16, default_value=1,
                           all_touched=True)

# Aggregate to 1 degree resolution
arr_1 = arr_0.reshape(arr_0.shape[0] // 10, 10,
                      arr_0.shape[1] // 10, 10).sum(axis=(1, 3))

print(arr_1)