python gdal完成arcgis分区统计功能（zonal）

最新推荐文章于 2024-05-06 01:30:29 发布
最新推荐文章于 2024-05-06 01:30:29 发布
阅读量3.4k 收藏 25
点赞数 4
分类专栏： gdal python
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_20373723/article/details/111350741
版权
python 同时被 2 个专栏收录
76 篇文章 24 订阅
订阅专栏
gdal
30 篇文章 30 订阅
订阅专栏

参考链接：https://towardsdatascience.com/zonal-statistics-algorithm-with-python-in-4-steps-382a3b66648a
注意事项，栅格需要和矢量的坐标系保持一致，结果存在了cvs文件中
结果
代码我稍微改了下，后续还会和矢量关联

import gdal
import ogr
import os
import numpy as np
import csv
import time

def boundingBoxToOffsets(bbox, geot):
    col1 = int((bbox[0] - geot[0]) / geot[1])
    col2 = int((bbox[1] - geot[0]) / geot[1]) + 1
    row1 = int((bbox[3] - geot[3]) / geot[5])
    row2 = int((bbox[2] - geot[3]) / geot[5]) + 1
    return [row1, row2, col1, col2]

def geotFromOffsets(row_offset, col_offset, geot):
    new_geot = [
    geot[0] + (col_offset * geot[1]),
    geot[1],
    0.0,
    geot[3] + (row_offset * geot[5]),
    0.0,
    geot[5]
    ]
    return new_geot

def setFeatureStats(fid, min, max, mean, median, sd, sum, count, names=["min", "max", "mean", "median", "sd", "sum", "count", "id"]):
    featstats = {
    names[0]: min,
    names[1]: max,
    names[2]: mean,
    names[3]: median,
    names[4]: sd,
    names[5]: sum,
    names[6]: count,
    names[7]: fid,
    }
    return featstats

def zonal(fn_raster, fn_zones, fn_csv):
    mem_driver = ogr.GetDriverByName("Memory")
    mem_driver_gdal = gdal.GetDriverByName("MEM")
    shp_name = "temp"

    # fn_raster = "C:/pyqgis/raster/USGS_NED_13_n45w116_IMG.img"
    # fn_zones = "C:/temp/zonal_stats/zones.shp"

    r_ds = gdal.Open(fn_raster)
    p_ds = ogr.Open(fn_zones)

    lyr = p_ds.GetLayer()
    geot = r_ds.GetGeoTransform()
    nodata = r_ds.GetRasterBand(1).GetNoDataValue()

    zstats = []

    p_feat = lyr.GetNextFeature()
    niter = 0

    while p_feat:
        if p_feat.GetGeometryRef() is not None:
            if os.path.exists(shp_name):
                mem_driver.DeleteDataSource(shp_name)
            tp_ds = mem_driver.CreateDataSource(shp_name)
            tp_lyr = tp_ds.CreateLayer('polygons', None, ogr.wkbPolygon)
            tp_lyr.CreateFeature(p_feat.Clone())
            offsets = boundingBoxToOffsets(p_feat.GetGeometryRef().GetEnvelope(),\
            geot)
            new_geot = geotFromOffsets(offsets[0], offsets[2], geot)

            tr_ds = mem_driver_gdal.Create(\
            "", \
            offsets[3] - offsets[2], \
            offsets[1] - offsets[0], \
            1, \
            gdal.GDT_Byte)

            tr_ds.SetGeoTransform(new_geot)
            gdal.RasterizeLayer(tr_ds, [1], tp_lyr, burn_values=[1])
            tr_array = tr_ds.ReadAsArray()

            r_array = r_ds.GetRasterBand(1).ReadAsArray(\
            offsets[2],\
            offsets[0],\
            offsets[3] - offsets[2],\
            offsets[1] - offsets[0])

            id = p_feat.GetFID()

            if r_array is not None:
                maskarray = np.ma.MaskedArray(\
                r_array,\
                maskarray=np.logical_or(r_array==nodata, np.logical_not(tr_array)))
                 
                if maskarray is not None:
                    zstats.append(setFeatureStats(\
                    id,\
                    maskarray.min(),\
                    maskarray.max(),\
                    maskarray.mean(),\
                    np.ma.median(maskarray),\
                    maskarray.std(),\
                    maskarray.sum(),\
                    maskarray.count()))
                else:
                    zstats.append(setFeatureStats(\
                    id,\
                    nodata,\
                    nodata,\
                    nodata,\
                    nodata,\
                    nodata,\
                    nodata,\
                    nodata))
            else:
                zstats.append(setFeatureStats(\
                    id,\
                    nodata,\
                    nodata,\
                    nodata,\
                    nodata,\
                    nodata,\
                    nodata,\
                    nodata))

            tp_ds = None
            tp_lyr = None
            tr_ds = None

            p_feat = lyr.GetNextFeature()

    # fn_csv = "C:/temp/zonal_stats/zstats.csv"
    col_names = zstats[0].keys()
    with open(fn_csv, 'w', newline='') as csvfile:
        writer = csv.DictWriter(csvfile, col_names)
        writer.writeheader()
        writer.writerows(zstats)

if __name__ == "__main__":
    time1 = time.time()
    fn_raster = './data/t_dem.tif'
    fn_zones = './data/grid1.shp'
    fn_csv = './data/zonsl.csv'
    zonal(fn_raster, fn_zones, fn_csv)
    time2 = time.time()
    print((time2-time1) / 3600.0)

如果只得到了csv而没有把值填进shp属性表那么上面操作的意义将大打折扣，下面是我完成以后的样子
对比图
统计结果
下面是更新后的代码：

import gdal
import ogr
import os
import numpy as np
import csv
import pandas as pd
import time

def boundingBoxToOffsets(bbox, geot):
    col1 = int((bbox[0] - geot[0]) / geot[1])
    col2 = int((bbox[1] - geot[0]) / geot[1]) + 1
    row1 = int((bbox[3] - geot[3]) / geot[5])
    row2 = int((bbox[2] - geot[3]) / geot[5]) + 1
    return [row1, row2, col1, col2]

def geotFromOffsets(row_offset, col_offset, geot):
    new_geot = [
    geot[0] + (col_offset * geot[1]),
    geot[1],
    0.0,
    geot[3] + (row_offset * geot[5]),
    0.0,
    geot[5]
    ]
    return new_geot

def setFeatureStats(fid, min, max, mean, median, sd, sum, count, names=["min", "max", "mean", "median", "sd", "sum", "count", "id"]):
    featstats = {
    names[0]: min,
    names[1]: max,
    names[2]: mean,
    names[3]: median,
    names[4]: sd,
    names[5]: sum,
    names[6]: count,
    names[7]: fid,
    }
    return featstats

def zonal(fn_raster, fn_zones, fn_csv):
    mem_driver = ogr.GetDriverByName("Memory")
    mem_driver_gdal = gdal.GetDriverByName("MEM")
    shp_name = "temp"

    # fn_raster = "C:/pyqgis/raster/USGS_NED_13_n45w116_IMG.img"
    # fn_zones = "C:/temp/zonal_stats/zones.shp"

    r_ds = gdal.Open(fn_raster)
    p_ds = ogr.Open(fn_zones)

    lyr = p_ds.GetLayer()
    geot = r_ds.GetGeoTransform()
    nodata = r_ds.GetRasterBand(1).GetNoDataValue()

    zstats = []

    p_feat = lyr.GetNextFeature()
    niter = 0

    while p_feat:
        if p_feat.GetGeometryRef() is not None:
            if os.path.exists(shp_name):
                mem_driver.DeleteDataSource(shp_name)
            tp_ds = mem_driver.CreateDataSource(shp_name)
            tp_lyr = tp_ds.CreateLayer('polygons', None, ogr.wkbPolygon)
            tp_lyr.CreateFeature(p_feat.Clone())
            offsets = boundingBoxToOffsets(p_feat.GetGeometryRef().GetEnvelope(),\
            geot)
            new_geot = geotFromOffsets(offsets[0], offsets[2], geot)

            tr_ds = mem_driver_gdal.Create(\
            "", \
            offsets[3] - offsets[2], \
            offsets[1] - offsets[0], \
            1, \
            gdal.GDT_Byte)

            tr_ds.SetGeoTransform(new_geot)
            gdal.RasterizeLayer(tr_ds, [1], tp_lyr, burn_values=[1])
            tr_array = tr_ds.ReadAsArray()

            r_array = r_ds.GetRasterBand(1).ReadAsArray(\
            offsets[2],\
            offsets[0],\
            offsets[3] - offsets[2],\
            offsets[1] - offsets[0])

            id = p_feat.GetFID()

            if r_array is not None:
                maskarray = np.ma.MaskedArray(\
                r_array,\
                maskarray=np.logical_or(r_array==nodata, np.logical_not(tr_array)))
                 
                if maskarray is not None:
                    zstats.append(setFeatureStats(\
                    id,\
                    maskarray.min(),\
                    maskarray.max(),\
                    maskarray.mean(),\
                    np.ma.median(maskarray),\
                    maskarray.std(),\
                    maskarray.sum(),\
                    maskarray.count()))
                else:
                    zstats.append(setFeatureStats(\
                    id,\
                    nodata,\
                    nodata,\
                    nodata,\
                    nodata,\
                    nodata,\
                    nodata,\
                    nodata))
            else:
                zstats.append(setFeatureStats(\
                    id,\
                    nodata,\
                    nodata,\
                    nodata,\
                    nodata,\
                    nodata,\
                    nodata,\
                    nodata))

            tp_ds = None
            tp_lyr = None
            tr_ds = None

            p_feat = lyr.GetNextFeature()
            
    col_names = zstats[0].keys()
    with open(fn_csv, 'w', newline='') as csvfile:
        writer = csv.DictWriter(csvfile, col_names)
        writer.writeheader()
        writer.writerows(zstats)

def shp_field_value(csv_file, shp):
    data = pd.DataFrame(pd.read_csv(csv_file))
    driver = ogr.GetDriverByName('ESRI Shapefile')
    layer_source = driver.Open(shp,1)
    lyr = layer_source.GetLayer()
    
    s_name = ogr.FieldDefn('min', ogr.OFTReal)
    lyr.CreateField(s_name)
    s_name = ogr.FieldDefn('max', ogr.OFTReal)
    lyr.CreateField(s_name)
    s_name = ogr.FieldDefn('mean', ogr.OFTReal)
    lyr.CreateField(s_name)
    s_name = ogr.FieldDefn('median', ogr.OFTReal)
    lyr.CreateField(s_name)
    s_name = ogr.FieldDefn('sd', ogr.OFTReal)
    lyr.CreateField(s_name)
    s_name = ogr.FieldDefn('sum', ogr.OFTReal)
    lyr.CreateField(s_name)
    s_name = ogr.FieldDefn('count', ogr.OFTReal)
    lyr.CreateField(s_name)

    count = 0
    defn = lyr.GetLayerDefn()
    featureCount = defn.GetFieldCount()
    feature = lyr.GetNextFeature()
    while feature is not None:
        for i in range(featureCount):
            feature.SetField('min', data['min'][count])
            feature.SetField('max', data['max'][count])
            feature.SetField('mean', data['mean'][count])
            feature.SetField('median', data['median'][count])
            feature.SetField('sd', data['sd'][count])
            feature.SetField('sum', data['sum'][count])
            feature.SetField('count', data['count'][count])
            lyr.SetFeature(feature)
        count+=1
        feature = lyr.GetNextFeature()

if __name__ == "__main__":
    time1 = time.time()
    fn_raster = './data/t_dem.tif'
    fn_zones = './data/grid1.shp'
    fn_csv = './data/zonsl.csv'
    zonal(fn_raster, fn_zones, fn_csv)
    shp_field_value(fn_csv, fn_zones)
    time2 = time.time()
    print((time2-time1) / 3600.0)

用时：0.00384633739789327 (h)

下面和arcgis的分区统计做下比较

import os
import time
import arcpy
from arcpy import env

def zonal(raster, shp):
	attri_table = "zonalstat"
	arcpy.gp.ZonalStatisticsAsTable_sa(shp, "FID", raster, attri_table, "NODATA", "ALL")
	arcpy.JoinField_management(shp,"FID", attri_table, "FID")

if __name__ == "__main__":
	time1 = time.time()
	dem_ras = './data/t_dem.tif'
	shp = './data/grid.shp'
	temp_path = './data/temp/'

	arcpy.CheckOutExtension('Spatial')
	arcpy.env.overwriteOutput = True
	env.workspace = temp_path

	zonal(dem_ras, shp)
	time2 = time.time()
	print((time2-time1) / 3600.0)

用时：0.196481666631 (h)
下面是arcgis统计以后的属性表，很明显arcgis多了几个字段包括AREA,MAJORITY,MINORITY等
arcgis
arcgis统计的字段比gdal实现的多，但是我感觉时间上还是有点太长了，就算gdal加上那些arcgis多出来的字段应该也不会这么长时间，感兴趣的朋友可以尝试一下。另外注意一下，两种方式得到的字段属性值也存在一定差异，但是大体上不会差很多，下面是mean值对应的视觉显示情况：
1.gdal
gdal

2.arcgis
arcgis
仔细观察下，可以发现gdal实现的效果还是有点瑕疵的，比如右下角的那几个异常值

测试数据链接：https://download.csdn.net/download/qq_20373723/13716488

如雾如电
关注
  • 4
    点赞
  • 踩
  • 25
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 8
    评论
专栏目录
gdal完成arcgis分区统计功能样例数据
12-18
gdal完成arcgis分区统计功能样例数据
python gdal安装与简单使用
09-18
Python GDAL库是一个强大的开源地理空间数据处理工具，它提供了对多种遥感和GIS数据格式的支持。GDAL（Geospatial Data Abstraction Library）最初是用C++编写的，但通过osgeo模块，我们可以使用Python接口来操作它...
ArcGIS应用（十九）Arcgis 统计分析计算多波段图像最大值、最小值、平均值等
soderayer的博客
05-29 8220
1.计算目的 计算多波段图像中，每个象元对应位置的所有波段值的最大值，最小值平均值等参数信息。 使用场景，比如计算时间序列遥感指数的最大值。假如有三个时间的NDVI指数，我们需要计算三个时间的最大值图像，就需要这样计算每个象元位置三个时相的最大值。 2.实现方法 2.1 使用工具 使用Arcgis统计分析工具计算得到，打开工具箱： 打开Spatial Analyst Tools中的zonal工具： 如下图： 使用ZonalStatistics工具进行统计分析.点击该工具，界面如下： .
【GDAL应用】基于rasterstats的矢量数据分区统计栅格值信息
最新发布
giser_rs
05-06 448
给定一个栅格数据和一个矢量面要素数据，计算矢量数据内栅格数据的统计值（如最大值、平均值、总和、最小值等）。返回一个列表，其中每个元素都是一个字典，表示与shp_path中相应多边形区域相关的统计信息。
ArcGIS教程：Zonal相关的工作原理及实例
qq_57970176的博客
11-02 4303
如果参数像元大小和环境像元大小均未指定，则默认的输出像元大小会根据输入数据集的类型，按照如下所示方式进行确定。如果区域输入和值输入均为具有相同像元大小的栅格，并且像元已对齐，则可将其直接应用于工具中，且工具执行期间不会对其进行内部重采样。如果区域输入和值输入均为具有相同像元大小的栅格，并且像元已对齐，则可将其直接应用于工具中，且工具执行期间不会对其进行内部重采样。如果区域输入和值输入均为具有相同像元大小的栅格，并且像元已对齐，则可将其直接应用于工具中，且工具执行期间不会对其进行内部重采样。
geemap学习笔记009：分区统计（zonal statistics）工具
主要是关于代码方面的内容，包括但不限于Python、MATLAB、ENVI、ArcGIS
11-20 635
分区统计（zonal statistics）工具
python 用geopandas，rasterio，gdal完成arcgis分区统计功能（zonal）
qq_20373723的博客
03-02 2699
前面已经有一篇实现分区统计功能的文章了（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是严格对应的，如果不对
gis里面python框_利用Python轻松在ArcGIS中进行图斑分类
weixin_29539581的博客
03-04 428
已知数据：图斑数据文件使用软件：ArcMap要求：国有土地与待确认土地，国有土地与待确认重叠部分超过待确认土地面积的一半则判定待确认土地为国有。下面开始操作：1、 打开ArcMap，将国有土地与待确认添加到数据框内，点击地理处理—相交输入要素待确认和国有土地，输出要素类位置选择我们自建的数据库，点击确定。2、 打开ArcToolbox—分析工具—统计分析—汇总统计数据输入表选择intersect，...
python分区统计批处理_python分区统计批处理_python_
10-03
6. **ArcGIS API for Python**：如果你的工作环境中使用了ArcGIS，那么可以利用其提供的Python API来进行分区统计。这个API包含了一个`ZonalStatisticsAsTable`函数，可以直接在Python脚本中调用，简化了处理流程。 ...
利用python GDAL库读写geotiff格式的遥感影像方法
09-19
Python中的GDAL库（Geospatial Data Abstraction Library）是处理这类数据的强大工具，它提供了读取、写入多种遥感和GIS文件格式的功能。本文将详细介绍如何使用Python GDAL库读写geotiff格式的遥感影像。 首先，要...
python3.9版本对应GDAL下载
07-27
python3.9版本对应GDAL下载
python GDAL实现按属性选择图层（等同于arcpy的SelectLayerByAttribute_management）
qq_20373723的博客
03-15 2271
用python版本的gdal完成按照字段属性选择矢量的功能
区域统计批处理功能（Python代码）
x2434417239的博客
12-18 1369
Tips ：本代码实现ArcGIS中的区域统计功能（Zonal），将统计的属性值直接赋予矢量文件属性表，可实现批量统计。 import ogr, os from rasterstats import zonal_stats def zonal(input_vector,in_raster,Stats_type,raster_name): shp = ogr.Open(input_vector,1) lyr = shp.GetLayer()
ArcGIS Zonal Statistics as Table 工具报错解决
半糖奶茶加珍珠
06-16 8240
该错误，主要是针对在使用Zonal Statistics as Table工具时修改了保存路径而发生的错误 在ArcGIS中执行Zonal Statistics as Table工具时，设置完参数后 （参数设置时，修改了存储位置，没有使用默认位置） 点击运行，出现如下图的错误： 解决方法： 1.点击File→Map Document Properties 2.在弹出的窗口中，看到最下面部分，在“Geodatabase”处，修改路径（点击那个红框中最右侧的文件夹图标） 3.该路径的选择，是你要执行工具时
基于Python的栅格数据分区统计
weixin_44785184的博客
04-01 1376
分区统计操作是一种用于计算由另一个数据集定义的区域内的栅格（值栅格）的像元值的统计操作。首先根据矢量数据将栅格化分区多个区域，提取每个区域的像元值并分别进行统计计算，而后输出结果（一般直接输出的输入矢量的某个字段里）。ArcGIS 有栅格的分区统计工具，但不够灵活，只能进行诸如最大值、最小值之类的统计计算。实现的栅格数据分区统计工具，能够更加灵活的实现分区统计功能，满足多样化的需求。，回复20230401获取示例数据和代码，包含这个工具的代码的写作思路。执行分区统计只需要调用类实例的。
python空间数据处理：part3-分区统计以及表格整理
weixin_45626690的博客
04-24 513
学习记录~~~
Python+GDAL实现矢量shape数据的分类统计
做GIS梦的人
12-29 664
之前在网上看到，并且也实践过，Python+GDAL可以实现对，主要是基于和。但是不知道如何实现分类统计的功能，查阅了很多资料都没有找到相关的API和案例。最后经过探索和实验，发现了基于，可以实现分类统计的功能，这里和大家分享一下。
利用GDAL对GRIB2数据进行矢栅叠加，并输出统计数据
做GIS梦的人
10-11 1166
遇到一个这样的需求。有一个包含多个天的GRIB2气象降雨数据，和一个河南省市级区划shp文件，需要统计每个市对每一天的最大降水值，并写入到Excel中，数据情况如下图所示：突然遇到这个问题，短时间内觉得还是比较麻烦的。
python gdal库
12-05
Python GDAL库是一个开源的地理数据抽象库。它提供了一种方便的方式来访问、读取和处理地理空间数据。GDAL库支持多种地理信息系统（GIS）格式，如Shapefile、GeoTIFF、KML等。 Python GDAL库的一个主要优势是它可以处理各种不同类型的地理数据并进行空间分析。它提供了强大的功能，如数据投影转换、裁剪、合并、重采样和地理空间分析等。 通过Python GDAL库，我们可以读取和写入地理矢量和栅格数据。例如，我们可以使用该库读取一个Shapefile文件，并将其转换为GeoJSON格式。我们还可以将一幅栅格图像裁剪为指定的区域，并保存为不同的格式。 Python GDAL库还可以进行地理空间分析。我们可以计算两个地理要素之间的距离，或者进行缓冲区分析，生成一定距离范围内的边界。此外，该库还支持地理要素之间的交叉、合并和裁剪等操作。 利用Python GDAL库，我们还可以进行地理数据的可视化。我们可以使用Matplotlib等可视化库将地理数据以图形的形式展示出来。这样可以更好地理解数据和展示结果。 总之，Python GDAL库是一个强大的工具，可用于读取、处理和分析各种地理空间数据。它提供了丰富的功能，同时易于使用，并且有大量的文档和示例代码可供参考。无论是进行地理数据处理、地理空间分析还是地理数据可视化，Python GDAL库都是一个不可或缺的工具。
评论 8
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

如雾如电

随缘

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值