在Python中用GDAL实现矢量对栅格的切割实例

概述：

本文讲述如何在Python中用GDAL实现根据输入矢量边界对栅格数据的裁剪。

效果：

裁剪前

矢量边界

裁剪后

实现代码：

# -*- coding: utf-8 -*-
"""
@author lzugis
@date 2017-06-02
@brief 利用shp裁剪影像
"""

from osgeo import gdal, gdalnumeric, ogr
from PIL import Image, ImageDraw
import os
import operator

gdal.UseExceptions()

# This function will convert the rasterized clipper shapefile
# to a mask for use within GDAL.
def imageToArray(i):
"""
Converts a Python Imaging Library array to a
gdalnumeric image.
"""
a=gdalnumeric.fromstring(i.tobytes(),'b')
a.shape=i.im.size[1], i.im.size[0]
return a

def arrayToImage(a):
"""
Converts a gdalnumeric array to a
Python Imaging Library Image.
"""
i=Image.frombytes('L',(a.shape[1],a.shape[0]),
(a.astype('b')).tobytes())
return i

def world2Pixel(geoMatrix, x, y):
"""
Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
the pixel location of a geospatial coordinate
"""
ulX = geoMatrix[0]
ulY = geoMatrix[3]
xDist = geoMatrix[1]
pixel = int((x - ulX) / xDist)
line = int((ulY - y) / xDist)
return (pixel, line)

#
# EDIT: this is basically an overloaded
# version of the gdal_array.OpenArray passing in xoff, yoff explicitly
# so we can pass these params off to CopyDatasetInfo
#
def OpenArray( array, prototype_ds = None, xoff=0, yoff=0 ):
ds = gdal.Open( gdalnumeric.GetArrayFilename(array) )

if ds is not None and prototype_ds is not None:
if type(prototype_ds).__name__ == 'str':
prototype_ds = gdal.Open( prototype_ds )
if prototype_ds is not None:
gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff )
return ds

def histogram(a, bins=range(0,256)):
"""
Histogram function for multi-dimensional array.
a = array
bins = range of numbers to match
"""
fa = a.flat
n = gdalnumeric.searchsorted(gdalnumeric.sort(fa), bins)
n = gdalnumeric.concatenate([n, [len(fa)]])
hist = n[1:]-n[:-1]
return hist

def stretch(a):
"""
Performs a histogram stretch on a gdalnumeric array image.
"""
hist = histogram(a)
im = arrayToImage(a)
lut = []
for b in range(0, len(hist), 256):
# step size
step = reduce(operator.add, hist[b:b+256]) / 255
# create equalization lookup table
n = 0
for i in range(256):
lut.append(n / step)
n = n + hist[i+b]
im = im.point(lut)
return imageToArray(im)

def main( shapefile_path, raster_path ):
# Load the source data as a gdalnumeric array
srcArray = gdalnumeric.LoadFile(raster_path)

# Also load as a gdal image to get geotransform
# (world file) info
srcImage = gdal.Open(raster_path)
geoTrans = srcImage.GetGeoTransform()

# Create an OGR layer from a boundary shapefile
shapef = ogr.Open(shapefile_path)
lyr = shapef.GetLayer( os.path.split( os.path.splitext( shapefile_path )[0] )[1] )
poly = lyr.GetNextFeature()

# Convert the layer extent to image pixel coordinates
minX, maxX, minY, maxY = lyr.GetExtent()
ulX, ulY = world2Pixel(geoTrans, minX, maxY)
lrX, lrY = world2Pixel(geoTrans, maxX, minY)

# Calculate the pixel size of the new image
pxWidth = int(lrX - ulX)
pxHeight = int(lrY - ulY)

clip = srcArray[:, ulY:lrY, ulX:lrX]

#
# EDIT: create pixel offset to pass to new image Projection info
#
xoffset = ulX
yoffset = ulY
print "Xoffset, Yoffset = ( %f, %f )" % ( xoffset, yoffset )

# Create a new geomatrix for the image
geoTrans = list(geoTrans)
geoTrans[0] = minX
geoTrans[3] = maxY

# Map points to pixels for drawing the
# boundary on a blank 8-bit,
# black and white, mask image.
points = []
pixels = []
geom = poly.GetGeometryRef()
pts = geom.GetGeometryRef(0)
for p in range(pts.GetPointCount()):
points.append((pts.GetX(p), pts.GetY(p)))
for p in points:
pixels.append(world2Pixel(geoTrans, p[0], p[1]))
rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)
rasterize = ImageDraw.Draw(rasterPoly)
rasterize.polygon(pixels, 0)
mask = imageToArray(rasterPoly)

# Clip the image using the mask
clip = gdalnumeric.choose(mask, \
(clip, 0)).astype(gdalnumeric.uint8)

# This image has 3 bands so we stretch each one to make them
# visually brighter
for i in range(3):
clip[i,:,:] = stretch(clip[i,:,:])

# Save new tiff
#
# EDIT: instead of SaveArray, let's break all the
# SaveArray steps out more explicity so
# we can overwrite the offset of the destination
# raster
#
### the old way using SaveArray
#
# gdalnumeric.SaveArray(clip, "OUTPUT.tif", format="GTiff", prototype=raster_path)
#
###
#
gtiffDriver = gdal.GetDriverByName( 'GTiff' )
if gtiffDriver is None:
raise ValueError("Can't find GeoTiff Driver")
gtiffDriver.CreateCopy( "beijing.tif",
OpenArray( clip, prototype_ds=raster_path, xoff=xoffset, yoff=yoffset )
)

# Save as an 8-bit jpeg for an easy, quick preview
clip = clip.astype(gdalnumeric.uint8)
gdalnumeric.SaveArray(clip, "beijing.jpg", format="JPEG")

gdal.ErrorReset()

if __name__ == '__main__':
#shapefile_path, raster_path
shapefile_path = 'beijing.shp'
raster_path = 'world.tif'
main( shapefile_path, raster_path )

补充知识：Python+GDAL | 读取矢量并写出txt

这篇文章主要描述了如何使用GDAL/OGR打开矢量文件、读取属性表，并将部分属性写出至txt。

代码

import ogr,sys,os
import numpy as np

os.chdir(r'E:\')

#设置driver,并打开矢量文件
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.Open('sites.shp', 0)
if ds is None:
print("Could not open", 'sites.shp')
sys.exit(1)
#获取图册
layer = ds.GetLayer()

#要素数量
numFeatures = layer.GetFeatureCount()
print("Feature count: "+str(numFeatures))

#获取范围
extent = layer.GetExtent()
print("Extent:", extent)
print("UL:", extent[0],extent[3])
print("LR:", extent[1],extent[2])

#获取要素
feature = layer.GetNextFeature()
ids = []
xs = []
ys = []
covers = []
#循环每个要素属性
while feature:
#获取字段“id”的属性
id = feature.GetField('id')
#获取空间属性
geometry = feature.GetGeometryRef()
x = geometry.GetX()
y = geometry.GetY()
cover = feature.GetField('cover')
ids.append(id)
xs.append(x)
ys.append(y)
covers.append(cover)
feature = layer.GetNextFeature()

data = [ids, xs, ys, covers]
data = np.array(data)
data = data.transpose()

#写出致txt
np.savetxt('myfile.txt',data, fmt='%s %s %s %s')
np.savetxt('myfile.csv',data, fmt='%s %s %s %s')

#释放文件空间
layer.ResetReading()
feature.Destroy()
ds.Destroy()

以上这篇在Python中用GDAL实现矢量对栅格的切割实例就是小编分享给大家的全部内容了，希望能给大家一个参考

您可能感兴趣的文章:

• GDAL 矢量属性数据修改方式(python)
• python+gdal+遥感图像拼接(mosaic)的实例
• python使用gdal对shp读取,新建和更新的实例