Rasterizing, separately

3draster.py

# PREPARATION
import geopandas as gp
from rasterio import features, Affine, float64
from rasterio.enums import MergeAlg
from xarray import DataArray
import numpy as np
import datashader as ds
from datashader import transfer_functions as tf, reductions as rd
from datashader.colors import Hot, viridis

import colorcet
kbc = colorcet.cm.kbc

def get_transform(cvs):
    '''transform into canvas coordinates'''
    return Affine((cvs.x_range[1] - cvs.x_range[0]) / cvs.plot_width, 0.0,
                   cvs.x_range[0], 0.0, (cvs.y_range[0] - cvs.y_range[1]) / cvs.plot_height, cvs.y_range[1])

cvs = ds.Canvas(plot_height=1000,
                plot_width=1000,
                x_range=(-73.647738, -74.301),
                y_range=(40.45, 41.))

aform = get_transform(cvs)

## ACTUAL PROCESS
gdf = pd.read_file(file.shp)


RS = []
def _rasterize_area(geom):
    raster = features.rasterize([(geom, geom.area)],  out_shape=(cvs.plot_height, cvs.plot_width), transform=aform, dtype=float64)
    RS.append(raster)
    
    
_ = gdf['geometry'].apply(_rasterize_area)


## NEXT STEPS
D = xarray.DataFRAME(RS)  # somehow convert list of 2d arrays into on2 3d-array
final = R.median(axis=3)
# and then visualize final