How to choose cell assignment method in GDAL/rasterio rasterize functions?

Often times, people (including myself) use gdal.rasterize() to distinguish data from no-data. It is a way to pluck a polygon or line which is in vector format out from there and into a raster. Essentially telling the raster where the alpha layer should be strong (transparent).

Now, you can still use this tool to do what you are wanting to do but if you want amazing control on the actual averaging of the precise data you have, you might want to consider using an array as an intermediate. Essentially, you take the vector data and push it into an array. Then you downgrade the size of the area to the corresponding granularity you want for the raster (1kmx1km per pixel) and then transform the array into a raster.

Since you are working with distinct values for the shapes, you don't want to mix these numbers together (fading green into pink) - you want hard boundaries. To do this, you can just write specific rules which propagate upwards (essentially mimicking the above method you want(which seems like you want the Maximum area method)). This is more work but it allows you to know exactly how the output is created and be able to back this up in writing.

I did take a quick look around and could not see anything concrete. Looking at this page from the GDAL FAQs, it looks like there used to be a dev page which explaining the tool in more detail but now this page throws a 404 (maybe it moved and you can find it?). However, right below that they do say

As Chris Barker suggests, GDAL's rasterization capability is not sophisticated from a render styling point of view. Other tools may be more appropriate if you want to do anything more sophisticated than rasterize the polygons in a single solid color.

And link to some possible other tools if you are interested.

I would say, if you are not afraid of working with arrays and writing simple rules to consolidate those arrays into a new size then take this route as you will have the most control

If I understand your core problem correctly you want to rasterize your shapefile and each pixel should be assign to the landcover class occupying the largest amount of area underneath each pixel.

I can present a solution which will not be the most elegant or cost effective solution processing wise (compared to the mentioned "array" solution), but with some "gridding" and multiprocessing (I will leave that part up to you) it might improve and its open source and it works.

It is a vector based solution mainly in geopandas and then rasterized in the end as a formality. I will explain some pseudo steps with code snippets.

Step 1: Create a vector grid over the extent of your LC shapefile (assuming you are working in a projected coordinate system with meters) in your required pixel size (e.g. 10 by 10m). Assign a GRID_ID column to the grid.

import numpy as np
from shapely.geometry import Polygon
import geopandas as gpd 

lc_shp_df = gpd.read_file('/path/to/your/LC/shape.shp') # with column e.g. 'class_code' for LC
xmin, ymin, xmax, ymax = lc_shp_df.unary_union.bounds

length= 10  # e.g. 10 by 10 meters
wide=   10
cols = list(range(int(np.floor(xmin)), int(np.ceil(xmax)), wide))
rows = list(range(int(np.floor(ymin+wide)), int(np.ceil(ymax+wide)), length))
rows.reverse()

polygons = []
for x in cols:
  for y in rows:
    polygons.append( Polygon([(x,y), (x+wide, y), (x+wide, y-length), (x, y-length)]) )
Grid_gdf = gpd.GeoDataFrame({'geometry':polygons})
Grid_gdf['GRID_ID'] = Grid_gdf.index

Step 2: Intersect the gird with the land cover shape and calculate the area of the intersection (note: ideally the land cover shape only has a landcover class column, this will make the process simpler, for the example I will use the column name class_code to represent the landcover class)

intersection = gpd.overlay(lc_shp_df, Grid_gdf, how='intersection')
intersection['area_msqr'] = intersection.area

In your created gird, for each GRID_ID you now have the area_msqr for each class_code

Step 3: Convert gpd to pd and get the class_code with the maximum area_msqr for each unique GRID_ID. Then join back to the original grid gpd and save as a shape.

import pandas as pd
intersection = pd.DataFrame(intersection)

# Find maximum class_code for each GIRD_ID
intersection_max = intersection.groupby('GRID_ID').idxmax()
class_code_max= intersection['class_code'].iloc[intersection_max['area_msqr']].reset_index()
class_code_max= class_code_max.rename(columns={'class_code':'max_cls_cd'})

# Join back to original gird dataframe
gird_crs = Grid_gdf.crs
Grid_gdf = pd.DataFrame(Grid_gdf)
Grid_gdf = Grid_gdf.set_index('GRID_ID')
LC_df = Grid_gdf.join(class_code_max['max_cls_cd'],how='left')
LC_df = gpd.GeoDataFrame(LC_df)
LC_df.crs = gird_crs 
LC_df.to_file('/path/to/your/Max_LC/shape.shp')

Step 4: Rasterize your new shape with the same geometrical properties (e.g. x,y boundaries, x y shape and 10 by 10 m) on the max_cls_cd attribute and you should have what you are looking for.

Since you have a 2GB LC vector file, I suggest to split this into parts and then batch execute. Hope that helps.