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
rasterstats. I would like to avoid using lower level libraries such as
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
calc_length_in_cell_listfunction 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.