栅格赋值矢量

最新推荐文章于 2024-06-12 17:45:28 发布

盐水鱼

最新推荐文章于 2024-06-12 17:45:28 发布

阅读量1.8k

点赞数 1

分类专栏： GDAL 文章标签： gdal python

本文链接：https://blog.csdn.net/weixin_43867016/article/details/107005865

版权

本文介绍了使用GDAL Python库进行栅格数据根据矢量数据赋值的方法，重点在于解决计算每个要素区域统计时可能出现的问题，如范围重叠导致的错误。通过按FID过滤图层来减少不必要的影响。

摘要由CSDN通过智能技术生成

栅格数据赋值矢量数据

前段时间给老师写栅格数据赋值矢量的代码，要用到GDAL包，网上搜了一遍没有找到合适的，只在GDAL的一个参考文档中找到一个类似的代码
Calculate zonal statistics

import gdal, ogr, osr, numpy
import sys


def zonal_stats(feat, input_zone_polygon, input_value_raster):

    # Open data
    raster = gdal.Open(input_value_raster)
    shp = ogr.Open(input_zone_polygon)
    lyr = shp.GetLayer()

    # Get raster georeference info
    transform = raster.GetGeoTransform()
    xOrigin = transform[0]
    yOrigin = transform[3]
    pixelWidth = transform[1]
    pixelHeight = transform[5]

    # Reproject vector geometry to same projection as raster
    sourceSR = lyr.GetSpatialRef()
    targetSR = osr.SpatialReference()
    targetSR.ImportFromWkt(raster.GetProjectionRef())
    coordTrans = osr.CoordinateTransformation(sourceSR,targetSR)
    feat = lyr.GetNextFeature()
    geom = feat.GetGeometryRef()
    geom.Transform(coordTrans)

    # Get extent of feat
    geom = feat.GetGeometryRef()
    if (geom.GetGeometryName() == 'MULTIPOLYGON'):
        count = 0
        pointsX = []; pointsY = []
        for polygon in geom:
            geomInner = geom.GetGeometryRef(count)
            ring = geomInner.GetGeometryRef(0)
            numpoints = ring.GetPointCount()
            for p in range(numpoints):
                    lon, lat, z = ring.GetPoint(p)
                    pointsX.append(lon)
                    pointsY.append(lat)
            count += 1
    elif (geom.GetGeometryName() == 'POLYGON'):
        ring = geom.GetGeometryRef(0)
        numpoints = ring