我编写了一个脚本来逐个光栅化不同的矢量文件，在此过程中，您必须输入该特定图层的每像素成本。最终，我试图实现一个包含所有初始矢量文件的光栅文件，这些文件现在都由单个光栅文件中的不同成本表示。
现在的问题是，当我输入的成本超过 255 时，光栅文件中的成本将仅为 255。这可能是因为 gdal 使用了 Byte 数据类型。我尝试将其更改为 GDAL_UInt16 以使 gat 值高达 65535，但这不起作用......当我使用以下代码输出文件时，我仍然会得到一个值在 0 到 255 之间的文件。
有谁能够帮我？
# -*- coding: utf-8 -*-
from osgeo import gdal, ogr, osr
import os
RASTERIZE_COST_FIELD = "__pixelkost__"
xmi = 0
xma = 0
ymi = 0
yma = 0
def rasterize(pixel_size=1):
# Ask for systempath towards main file
print ""
print "**************************************************"
print "** **"
print "** Please enter the systempath towards the **"
print "** shapefile which ranges over the total sur- **"
print "** face area of the used region. This will **"
print "** probably be the file with administrative **"
print "** lots. **"
print "** **"
print "**************************************************"
print ""
path = one()
# Ask for cost
cost1 = three()
# produce the outputpath
(shapefilefilepath, shapefilename) = os.path.split(path)
(shapefileshortname, extension) = os.path.splitext(shapefilename)
ofp = shapefilefilepath + "/" + shapefileshortname + ".img"
# Open the data source
orig_data_source = ogr.Open(path)
# Make a copy of the layer's data source because we'll need to
# modify its attributes table
source_ds = ogr.GetDriverByName("Memory").CopyDataSource(
orig_data_source, "")
source_layer = source_ds.GetLayer(0)
source_srs = source_layer.GetSpatialRef()
x_min, x_max, y_min, y_max = source_layer.GetExtent()
xmi = x_min
ymi = y_min
xma = x_max
yma = y_max
# Create a field in the source layer to hold the features pixelkost
field_def = ogr.FieldDefn(RASTERIZE_COST_FIELD, ogr.OFTReal)
source_layer.CreateField(field_def)
source_layer_def = source_layer.GetLayerDefn()
field_index = source_layer_def.GetFieldIndex(RASTERIZE_COST_FIELD)
# Generate random values for the field (it's here that the value
# of the attribute should be used, but you get the idea)
for feature in source_layer:
feature.SetField(field_index, cost1)
source_layer.SetFeature(feature)
# Create the destination data source
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('GTiff').Create(ofp, x_res,
y_res, 1, gdal.GDT_UInt16)
target_ds.SetGeoTransform((
x_min, pixel_size, 0,
y_max, 0, -pixel_size,
))
if source_srs:
# Make the target raster have the same projection as the source
target_ds.SetProjection(source_srs.ExportToWkt())
else:
# Source has no projection (needs GDAL >= 1.7.0 to work)
target_ds.SetProjection(ofp)
# Rasterize
err = gdal.RasterizeLayer(target_ds, [1], source_layer,
burn_values=[2509],
options=["ATTRIBUTE=%s" % RASTERIZE_COST_FIELD])
target_ds = None
Array1 = raster2array(ofp)
a = 1
print ""
print "**************************************************"
print "** **"
print "** Please enter, one by one, the systempath **"
print "** towards the shapefile you want to add to **"
print "** the set of rasterfiles. These will be **"
print "** merged to one final cost-file. **"
print "** **"
print "** If you want to stop adding shapefiles, type **"
print "** STOP **"
print "** **"
print "**************************************************"
print ""
while a != "STOP":
# Ask for systempath towards next file
path2 = one()
# Ask for cost
cost2 = three()
print ""
print "This process may take a while, please be patient."
# produce the outputpath
(shapefilefilepath, shapefilename) = os.path.split(path2)
(shapefileshortname, extension) = os.path.splitext(shapefilename)
ofp2 = shapefilefilepath + "/" + shapefileshortname + ".img"
# Open the data source
orig_data_source = ogr.Open(path2)
# Make a copy of the layer's data source because we'll need to
# modify its attributes table
source_ds = ogr.GetDriverByName("Memory").CopyDataSource(
orig_data_source, "")
source_layer = source_ds.GetLayer(0)
source_srs = source_layer.GetSpatialRef()
# Create a field in the source layer to hold the features pixelkost
field_def = ogr.FieldDefn(RASTERIZE_COST_FIELD, ogr.OFTReal)
source_layer.CreateField(field_def)
source_layer_def = source_layer.GetLayerDefn()
field_index = source_layer_def.GetFieldIndex(RASTERIZE_COST_FIELD)
for feature in source_layer:
feature.SetField(field_index, cost2)
source_layer.SetFeature(feature)
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('GTiff').Create(ofp2, x_res,
y_res, 1, gdal.GDT_UInt16)
target_ds.SetGeoTransform((
x_min, pixel_size, 0,
y_max, 0, -pixel_size,
))
if source_srs:
target_ds.SetProjection(source_srs.ExportToWkt())
else:
target_ds.SetProjection(ofp2)
err = gdal.RasterizeLayer(target_ds, [1], source_layer,
burn_values=[2509],
options=["ATTRIBUTE=%s" % RASTERIZE_COST_FIELD])
target_ds = None
Array2 = raster2array(ofp2)
length = len(Array1)
length2 = len(Array1[1])
l = 0
while l < length:
l2 = 0
while l2 < length2:
Array1[l,l2] += Array2[l,l2]
l2 += 1
l2 = 0
l += 1
while True:
try:
a = str(input("Press enter if you want to continue inputting shapefiles. Otherwise type STOP: "))
break
except (NameError, SyntaxError, ValueError, TypeError):
print ""
print("This is not a string. Please do not forget to type your input between quotation marks.")
print ""
except (UnicodeEncodeError):
print("Please do not use special characters in the systempath.")
print ""
return (Array1, path, ofp)
def one():
while True:
try:
syspath = str(input("Systempath: "))
break
except (NameError, SyntaxError, ValueError, TypeError):
print ""
print("This is not a string. Please do not forget to type your input between quotation marks.")
print ""
except (UnicodeEncodeError):
print("Please do not use special characters in the systempath.")
print ""
return syspath
def three():
while True:
try:
cost = int(input("Cost per squared meter for this type of land-use: "))
break
except (NameError, ValueError, SyntaxError, UnicodeEncodeError):
print ""
print("This is not a numerical value, please re-enter.")
print ""
return cost
def raster2array(rasterfn):
raster = gdal.Open(rasterfn)
band = raster.GetRasterBand(1)
array = band.ReadAsArray()
return array
def array2raster(newRasterfn,rasterfn,array):
raster = gdal.Open(rasterfn)
geotransform = raster.GetGeoTransform()
originX = geotransform[0] # East/West location of Upper Left corner
originY = geotransform[3] # North/South location of Upper Left corner
pixelWidth = geotransform[1] # X pixel size
pixelHeight = geotransform[5] # Y pixel size
cols = array.shape[1]
rows = array.shape[0]
driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create(newRasterfn, cols, rows, gdal.GDT_UInt16)
outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
outband = outRaster.GetRasterBand(1)
outband.WriteArray(array)
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
outRaster.SetProjection(outRasterSRS.ExportToWkt())
outband.FlushCache()
if __name__ == '__main__':
A1, pa, imgpa = rasterize()
(shapefilefilepath, shapefilename) = os.path.split(pa)
ofp3 = shapefilefilepath + "/merge.img"
array2raster(ofp3,imgpa,A1)
print ""
print "**************************************************"
print ""
print "The merged rasterfile has been saved at:"
print ""
print str(ofp3)
print ""
print "**************************************************"