前面已经有一篇实现分区统计功能的文章了(https://blog.csdn.net/qq_20373723/article/details/111350741),但是核心计算部分是用的numpy,比较难实现majority这个属性,这里用了geopandas,rasterio,gdal地理处理模块组合实现了自动分区统计到表,希望对大家有所帮助
python 3.7.3,geopandas 0.8.2 ,GDAL 3.2.1,rasterio 1.2.0,其中的gdal和rasterio是严格对应的,如果不对应,rasterio是无法安装的。
结果为证:
代码:
import geopandas as pd
import rasterio
from rasterio.plot import show
from rasterstats import zonal_stats
from osgeo import gdal, ogr
ras_path = './fixed_0.tif'
shp_path = './v2.shp'
stats_list = ['min', 'max', 'mean', 'median', '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)
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()
参考链接:https://pythonhosted.org/rasterstats/manual.html#vector-data-sources
上面的代码可能会遇到异常,第一次运行属性表没有改变,下面的代码应该不会有问题
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)