在Python中用GDAL实现矢量对栅格的切割实例
概述：
本文讲述如何在Python中用GDAL实现根据输入矢量边界对栅格数据的裁剪。
效果：
裁剪前
矢量边界
裁剪后
实现代码：
#-*-coding:utf-8-*-
"""
@authorlzugis
@date2017-06-02
@brief利用shp裁剪影像
"""
fromosgeoimportgdal,gdalnumeric,ogr
fromPILimportImage,ImageDraw
importos
importoperator
gdal.UseExceptions()
#Thisfunctionwillconverttherasterizedclippershapefile
#toamaskforusewithinGDAL.
defimageToArray(i):
"""
ConvertsaPythonImagingLibraryarraytoa
gdalnumericimage.
"""
a=gdalnumeric.fromstring(i.tobytes(),'b')
a.shape=i.im.size[1],i.im.size[0]
returna
defarrayToImage(a):
"""
Convertsagdalnumericarraytoa
PythonImagingLibraryImage.
"""
i=Image.frombytes('L',(a.shape[1],a.shape[0]),
(a.astype('b')).tobytes())
returni
defworld2Pixel(geoMatrix,x,y):
"""
Usesagdalgeomatrix(gdal.GetGeoTransform())tocalculate
thepixellocationofageospatialcoordinate
"""
ulX=geoMatrix[0]
ulY=geoMatrix[3]
xDist=geoMatrix[1]
pixel=int((x-ulX)/xDist)
line=int((ulY-y)/xDist)
return(pixel,line)
#
#EDIT:thisisbasicallyanoverloaded
#versionofthegdal_array.OpenArraypassinginxoff,yoffexplicitly
#sowecanpasstheseparamsofftoCopyDatasetInfo
#
defOpenArray(array,prototype_ds=None,xoff=0,yoff=0):
ds=gdal.Open(gdalnumeric.GetArrayFilename(array))
ifdsisnotNoneandprototype_dsisnotNone:
iftype(prototype_ds).__name__=='str':
prototype_ds=gdal.Open(prototype_ds)
ifprototype_dsisnotNone:
gdalnumeric.CopyDatasetInfo(prototype_ds,ds,xoff=xoff,yoff=yoff)
returnds
defhistogram(a,bins=range(0,256)):
"""
Histogramfunctionformulti-dimensionalarray.
a=array
bins=rangeofnumberstomatch
"""
fa=a.flat
n=gdalnumeric.searchsorted(gdalnumeric.sort(fa),bins)
n=gdalnumeric.concatenate([n,[len(fa)]])
hist=n[1:]-n[:-1]
returnhist
defstretch(a):
"""
Performsahistogramstretchonagdalnumericarrayimage.
"""
hist=histogram(a)
im=arrayToImage(a)
lut=[]
forbinrange(0,len(hist),256):
#stepsize
step=reduce(operator.add,hist[b:b+256])/255
#createequalizationlookuptable
n=0
foriinrange(256):
lut.append(n/step)
n=n+hist[i+b]
im=im.point(lut)
returnimageToArray(im)
defmain(shapefile_path,raster_path):
#Loadthesourcedataasagdalnumericarray
srcArray=gdalnumeric.LoadFile(raster_path)
#Alsoloadasagdalimagetogetgeotransform
#(worldfile)info
srcImage=gdal.Open(raster_path)
geoTrans=srcImage.GetGeoTransform()
#CreateanOGRlayerfromaboundaryshapefile
shapef=ogr.Open(shapefile_path)
lyr=shapef.GetLayer(os.path.split(os.path.splitext(shapefile_path)[0])[1])
poly=lyr.GetNextFeature()
#Convertthelayerextenttoimagepixelcoordinates
minX,maxX,minY,maxY=lyr.GetExtent()
ulX,ulY=world2Pixel(geoTrans,minX,maxY)
lrX,lrY=world2Pixel(geoTrans,maxX,minY)
#Calculatethepixelsizeofthenewimage
pxWidth=int(lrX-ulX)
pxHeight=int(lrY-ulY)
clip=srcArray[:,ulY:lrY,ulX:lrX]
#
#EDIT:createpixeloffsettopasstonewimageProjectioninfo
#
xoffset=ulX
yoffset=ulY
print"Xoffset,Yoffset=(%f,%f)"%(xoffset,yoffset)
#Createanewgeomatrixfortheimage
geoTrans=list(geoTrans)
geoTrans[0]=minX
geoTrans[3]=maxY
#Mappointstopixelsfordrawingthe
#boundaryonablank8-bit,
#blackandwhite,maskimage.
points=[]
pixels=[]
geom=poly.GetGeometryRef()
pts=geom.GetGeometryRef(0)
forpinrange(pts.GetPointCount()):
points.append((pts.GetX(p),pts.GetY(p)))
forpinpoints:
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)
#Cliptheimageusingthemask
clip=gdalnumeric.choose(mask,\
(clip,0)).astype(gdalnumeric.uint8)
#Thisimagehas3bandssowestretcheachonetomakethem
#visuallybrighter
foriinrange(3):
clip[i,:,:]=stretch(clip[i,:,:])
#Savenewtiff
#
#EDIT:insteadofSaveArray,let'sbreakallthe
#SaveArraystepsoutmoreexplicityso
#wecanoverwritetheoffsetofthedestination
#raster
#
###theoldwayusingSaveArray
#
#gdalnumeric.SaveArray(clip,"OUTPUT.tif",format="GTiff",prototype=raster_path)
#
###
#
gtiffDriver=gdal.GetDriverByName('GTiff')
ifgtiffDriverisNone:
raiseValueError("Can'tfindGeoTiffDriver")
gtiffDriver.CreateCopy("beijing.tif",
OpenArray(clip,prototype_ds=raster_path,xoff=xoffset,yoff=yoffset)
)
#Saveasan8-bitjpegforaneasy,quickpreview
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。
代码
importogr,sys,os
importnumpyasnp
os.chdir(r'E:\')
#设置driver,并打开矢量文件
driver=ogr.GetDriverByName('ESRIShapefile')
ds=driver.Open('sites.shp',0)
ifdsisNone:
print("Couldnotopen",'sites.shp')
sys.exit(1)
#获取图册
layer=ds.GetLayer()
#要素数量
numFeatures=layer.GetFeatureCount()
print("Featurecount:"+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=[]
#循环每个要素属性
whilefeature:
#获取字段“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实现矢量对栅格的切割实例就是小编分享给大家的全部内容了，希望能给大家一个参考，也希望大家多多支持毛票票。
声明：本文内容来源于网络，版权归原作者所有，内容由互联网用户自发贡献自行上传，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任。如果您发现有涉嫌版权的内容，欢迎发送邮件至：czq8825#qq.com（发邮件时，请将#更换为@）进行举报，并提供相关证据，一经查实，本站将立刻删除涉嫌侵权内容。