Как можно преобразовать линейные векторы в растр линейных плотностей?

Question

Я пытаюсь преобразовать векторный слой из нескольких строк в растр, где значения пикселей равны количеству линий, проходящих через эту ячейку.

Я не уверен, с чего начать. Любая помощь будет оценена. Спасибо!

Answer 1

import geopandas as gpd
import numpy as np
import pandas as pd
from shapely.geometry import Polygon
from shapely.geometry import LineString
from geopandas import sjoin
import rasterio
import xarray as xr
from geocube.api.core import make_geocube
import os

# make some linestrings
line1 = LineString([(17.45, 51.11), (17.46, 51.35)])
line2 = LineString([(17.2, 51.8), (17.23, 51.72)])
line3 = LineString([(16.91, 51.3), (16.84, 50.65)])
line4 = LineString([(16.31, 51.44), (16.37, 51.4)])
line5 = LineString([(16.38, 51.42), (16.43, 51.49)])

# make linestring gdf and save to .shp
line_gdf = gpd.GeoDataFrame([line1, line2, line3, line4, line5],
                            geometry=[line1, line2, line3, line4, line5])

line_gdf.crs = "EPSG:4326"
line_gdf = line_gdf[['geometry']]
line_gdf.to_file('line_gdf.shp')

# make shapely polygon grid 
xmin,ymin,xmax,ymax = [16, 50, 18, 52]

length = 0.25
width = 0.25

cols = list(np.arange(np.floor(xmin), np.ceil(xmax), width))
rows = list(np.arange(np.floor(ymin), np.ceil(ymax), length))
rows.reverse()

polygons = []
for x in cols:
    for y in rows:
        polygons.append( Polygon([(x,y), (x+width, y), (x+width, y-length), (x, y-length)]) )

#make grid_gdf, save to .shp
grid_gdf = gpd.GeoDataFrame({'geometry':polygons})
grid_gdf.crs = "EPSG:4326"
grid_gdf.to_file('grid_gdf.shp')

# loop through grid cells and count how many lines intersect each grid cell
int_count_dict = {}
for idx, row in grid_gdf.iterrows():
    lines_in_polygons_gdf = gpd.sjoin(line_gdf,
                                      grid_gdf.iloc[idx:idx+1,:],
                                      op='intersects',
                                      how='inner')
    # store line intersections as dict values
    int_count_dict[idx] = lines_in_polygons_gdf.shape[0]

#map intersections to grid_gdf into count column
grid_gdf['count'] = grid_gdf.index.map(int_count_dict)

#make xarray of grid_gdf
cube = make_geocube(vector_data=grid_gdf, measurements=['count'], resolution=(0.25, -0.25))

#export to .tif
cube["count"].rio.to_raster("line_count_raster.tif")

print('gpd version:', gpd.__version__)
print('np version:', np.__version__)
print('pd version:', pd.__version__)
print('shapely version:', shapely.__version__)
print('rasterio version:', rasterio.__version__)
print('xarray version:', xr.__version__)
print('geocube version:', geocube.__version__)

gpd version: 0.7.0
np version: 1.18.1
pd version: 1.0.3
shapely version: 1.7.0
rasterio version: 1.1.3
xarray version: 0.15.1
geocube version: 0.0.11