概要

本文主要利用了majority这个字段，来区别栅格与矢量ploygon_feature之间的关系。简单来说，当ploygon_feature内对应的栅格区域，栅格的某个值X超过了百分之五十，那么majory这个字段对应的值就为X。

代码

import time
import geopandas as pd
import rasterio
from rasterio.plot import show
from rasterstats import zonal_stats
from osgeo import gdal, ogr

start = time.time()
ras_path = 'cq_test.tif'
shp_path = 'cq_test.shp'
# stats_list = ['min', 'max', 'mean', 'median', 'majority']
stats_list = ['majority']

ras_driver = rasterio.open(ras_path)
array = ras_driver.read(1)
affine = ras_driver.transform
shp_driver = pd.read_file(shp_path)
zs = zonal_stats(shp_path, array, affine=affine, stats = stats_list)

driver = ogr.GetDriverByName('ESRI Shapefile')
layer_source = driver.Open(shp_path,1)
lyr = layer_source.GetLayer()
defn = lyr.GetLayerDefn()

featureCount = defn.GetFieldCount()
exists_fields = []
for i in range(featureCount):
	field = defn.GetFieldDefn(i)
	field_name = field.GetNameRef()
	exists_fields.append(field_name)

for ele in stats_list:
	if ele in exists_fields:
		pass
	else:
		# cls_name = ogr.FieldDefn(k, ogr.OFTString)
		cls_name = ogr.FieldDefn(ele, ogr.OFTReal)
		# cls_name.SetWidth(64)
		lyr.CreateField(cls_name)

driver = None

driver = ogr.GetDriverByName('ESRI Shapefile')
layer_source = driver.Open(shp_path,1)
lyr = layer_source.GetLayer()
defn = lyr.GetLayerDefn()

featureCount = defn.GetFieldCount()

count = 0
feature = lyr.GetNextFeature()
while feature is not None:
	for i in range(featureCount):
		field = defn.GetFieldDefn(i)
		field_name = field.GetNameRef()
		if field_name in stats_list:
			feature.SetField(field_name, zs[count][field_name])
			lyr.SetFeature(feature)
		else:
			pass
	count+=1
	feature = lyr.GetNextFeature()

end = time.time()
print((end-start)/3600.0)