为什么使用 GDAL 在 Python 中进行光栅化不起作用?
Why rasterizing in Python using GDAL does not work?
我正在使用 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)
我使用 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 中读取包含 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)
我使用 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 时,不要以米为单位指定...