0

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:

enter image description here

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 enter image description here 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.

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)
  • Thanks for your work. This method is not unfamiliar to me. However it does not measure the lengths accurately enough for my project, furthermore if forces you to aggregate the "finer resolution" to something which is an exact multiple of the "coarser resolution". That is to say, the result output by your script will be accurate to 0.1 degrees but not more. Unfortunately this will not cut it for our particular project.... – user32882 Oct 11 '18 at 18:37

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.