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

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    

并始终记住保持像元大小与坐标系相关，例如不要'当 shapefile 的投影为 WGS 时，t 以米为单位指定...