为什么使用 GDAL 在 Python 中进行光栅化不起作用？

Question

我正在使用 GDAL 在 Python 中读取包含 0 到 100 范围内数据的 shapefile。不幸的是，虽然它没有给出错误，但结果是不正确的（与 QGIS 相比）。我尝试了不同的 NoDataValue，但没有找到正确的结果。

代码如下：

from osgeo import gdal
from osgeo import ogr
import matplotlib.pyplot as plt
import numpy as np
import glob
import numpy.ma as ma

def Feature_to_Raster(input_shp, output_tiff, cellsize, field_name=True, NoData_value=-9999):
    # Input
    inp_driver = ogr.GetDriverByName('ESRI Shapefile')
    inp_source = inp_driver.Open(input_shp, 0)
    inp_lyr = inp_source.GetLayer(0)
    inp_srs = inp_lyr.GetSpatialRef()

    # Extent
    x_min, x_max, y_min, y_max = inp_lyr.GetExtent()
    x_ncells = int((x_max - x_min) / cellsize)
    y_ncells = int((y_max - y_min) / cellsize)

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,1, gdal.GDT_Float32)
    out_source.SetGeoTransform((x_min, cellsize, 0, y_max, 0, -cellsize))
    out_source.SetProjection(inp_srs.ExportToWkt())
    out_lyr = out_source.GetRasterBand(1)
    out_lyr.SetNoDataValue(NoData_value)

    # Rasterize
    # print(inp_lyr)
    if field_name:
        gdal.RasterizeLayer(out_source, [1], inp_lyr, options=["ATTRIBUTE=CT"])
    else:
        gdal.RasterizeLayer(out_source, [1], inp_lyr, burn_values=[1])
    
   
     # Save and/or close the data sources
    inp_source = None
    out_source = None
 
    
    ds= gdal.Open('name.tif')
    ndv= ds.GetRasterBand(1).GetNoDataValue()
    bnd1= ds.GetRasterBand(1).ReadAsArray()
    bnd1[bnd1==ndv]= np.nan
    tt= ma.masked_outside(bnd1, 1,100)
    plt.imshow(tt, cmap='jet')
    plt.colorbar()
    plt.xlabel('Column #')
    plt.ylabel('Row #')
    plt.show()    

    # Return
    return output_tiff
    
output_tiff= 'D:/myfolder/name.tif'
input_shp= 'D:/myfolder/cis_SGRDAMID_20101201.shp'
Feature_to_Raster(input_shp, output_tiff, cellsize, field_name=True, NoData_value=-9999)

Answer 1

我使用 gdal.Rasterize 函数取得了更大的成功

看看这是否能解决您的问题：

你可以替换这个：

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,1, gdal.GDT_Float32)
    out_source.SetGeoTransform((x_min, cellsize, 0, y_max, 0, -cellsize))
    out_source.SetProjection(inp_srs.ExportToWkt())
    out_lyr = out_source.GetRasterBand(1)
    out_lyr.SetNoDataValue(NoData_value)

    # Rasterize
    # print(inp_lyr)
    if field_name:
        gdal.RasterizeLayer(out_source, [1], inp_lyr, options=["ATTRIBUTE=CT"])
    else:
        gdal.RasterizeLayer(out_source, [1], inp_lyr, burn_values=[1])

有了这个：

if field_name:
    # This will rasterize your shape file according to the specified attribute field
    rasDs = gdal.Rasterize(output_tiff, input_shp,
               xRes=cellsize, yRes=cellsize,
               outputBounds=[x_min, y_min,x_max, y_max],
               noData=NoData_value,
               outputType=gdal.GDT_Float32
               attribute='CT', # Or whatever your attribute field name is
               allTouched=True)
else:
    # This will just give 255 where there are vector data since no attribute is defined
    rasDs = gdal.Rasterize(output_tiff, input_shp,
               xRes=cellsize, yRes=cellsize,
               outputBounds=[x_min, y_min,x_max, y_max],
               noData=NoData_value,
               outputType=gdal.GDT_Float32
               allTouched=True)

rasDs = inp_source = None

并始终记得保持您的 cell-size 与您的坐标系相关，例如当 shapefile 的投影是 WGS 时，不要以米为单位指定...

为什么使用 GDAL 在 Python 中进行光栅化不起作用？

Why rasterizing in Python using GDAL does not work?

python

gdal

shapefile

qgis

rasterizing