利用Python裁切tiff图像且读取tiff,shp文件的实例
我就废话不多说了，还是直接看代码吧！
fromosgeoimportgdal,gdalnumeric,ogr
fromPILimportImage,ImageDraw
fromosgeoimportgdal_array
importos
importoperator
fromfunctoolsimportreduce
gdal.UseExceptions()
defreadTif(fileName):
dataset=gdal.Open(fileName)
ifdataset==None:
print(fileName+"文件无法打开")
return
im_width=dataset.RasterXSize#栅格矩阵的列数
im_height=dataset.RasterYSize#栅格矩阵的行数
im_bands=dataset.RasterCount#波段数
band1=dataset.GetRasterBand(1)
print(band1)
print('BandType=',gdal.GetDataTypeName(band1.DataType))
im_data=dataset.ReadAsArray(0,0,im_width,im_height)#获取数据
im_geotrans=dataset.GetGeoTransform()#获取仿射矩阵信息
im_proj=dataset.GetProjection()#获取投影信息
im_blueBand=im_data[0,0:im_height,0:im_width]#获取蓝波段
im_greenBand=im_data[1,0:im_height,0:im_width]#获取绿波段
im_redBand=im_data[2,0:im_height,0:im_width]#获取红波段
im_nirBand=im_data[3,0:im_height,0:im_width]#获取近红外波段
return(im_width,im_height,im_bands,im_data,im_geotrans
,im_proj,im_blueBand,im_greenBand,im_redBand,im_nirBand)
#保存tif文件函数
importgdal
importnumpyasnp
defwriteTiff(im_data,im_width,im_height,im_bands,im_geotrans,im_proj,path):
if'int8'inim_data.dtype.name:
datatype=gdal.GDT_Byte
elif'int16'inim_data.dtype.name:
datatype=gdal.GDT_UInt16
else:
datatype=gdal.GDT_Float32
iflen(im_data.shape)==3:
im_bands,im_height,im_width=im_data.shape
eliflen(im_data.shape)==2:
im_data=np.array([im_data])
else:
im_bands,(im_height,im_width)=1,im_data.shape
#创建文件
driver=gdal.GetDriverByName("GTiff")
dataset=driver.Create(path,im_width,im_height,im_bands,datatype)
if(dataset!=None):
dataset.SetGeoTransform(im_geotrans)#写入仿射变换参数
dataset.SetProjection(im_proj)#写入投影
foriinrange(im_bands):
dataset.GetRasterBand(i+1).WriteArray(im_data[i])
deldataset
#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_array.OpenArray(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(4):
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("beijing1.tif",
OpenArray(clip,prototype_ds=raster_path,xoff=xoffset,yoff=yoffset)
)
print(raster_path)
#Saveasan8-bitjpegforaneasy,quickpreview
clip=clip.astype(gdalnumeric.uint8)
gdalnumeric.SaveArray(clip,"beijing1.jpg",format="JPEG")
gdal.ErrorReset()
if__name__=='__main__':
#shapefile_path,raster_path
shapefile_path=r'C:\Users\Administrator\Desktop\裁切shp\New_Shapefile.shp'
raster_path=r'C:\Users\Administrator\Desktop\2230542.tiff'
main(shapefile_path,raster_path)
补充知识：python代码裁剪tiff影像图和转换成png格式+裁剪Png图片
先来看一下需要转换的tiff原始图的信息，如下图所示。
tiff转换成png和裁剪tiff的代码（opencv）
importcv2ascv importos """ 转换tiff格式为png+横向裁剪tiff遥感影像图 """ defConvert_To_Png_AndCut(dir): files=os.listdir(dir) ResultPath1="./RS_ToPngDir/"#定义转换格式后的保存路径 ResultPath2="./RS_Cut_Result/"#定义裁剪后的保存路径 ResultPath3="./RS_Cut_Result/"#定义裁剪后的保存路径 forfileinfiles:#这里可以去掉for循环 a,b=os.path.splitext(file)#拆分影像图的文件名称 this_dir=os.path.join(dir+file)#构建保存路径+文件名 img=cv.imread(this_dir,1)#读取tif影像 #第二个参数是通道数和位深的参数， #IMREAD_UNCHANGED=-1#不进行转化，比如保存为了16位的图片，读取出来仍然为16位。 #IMREAD_GRAYSCALE=0#进行转化为灰度图，比如保存为了16位的图片，读取出来为8位，类型为CV_8UC1。 #IMREAD_COLOR=1#进行转化为RGB三通道图像，图像深度转为8位 #IMREAD_ANYDEPTH=2#保持图像深度不变，进行转化为灰度图。 #IMREAD_ANYCOLOR=4#若图像通道数小于等于3，则保持原通道数不变；若通道数大于3则只取取前三个通道。图像深度转为8位 cv.imwrite(ResultPath1+a+"_"+".png",img)#保存为png格式 #下面开始裁剪-不需要裁剪tiff格式的可以直接注释掉 hight=img.shape[0]#opencv写法，获取宽和高 width=img.shape[1] #定义裁剪尺寸 w=480#宽度 h=360#高度 _id=1#裁剪结果保存文件名：0-N升序方式 i=0 while(i+h<=hight):#控制高度,图像多余固定尺寸总和部分不要了 j=0 while(j+w<=width):#控制宽度，图像多余固定尺寸总和部分不要了 cropped=img[i:i+h,j:j+w]#裁剪坐标为[y0:y1,x0:x1] cv.imwrite(ResultPath2+a+"_"+str(_id)+b,cropped) _id+=1 j+=w i=i+h """ 横向裁剪PNG图 """ deftoCutPng(dir): files=os.listdir(dir) ResultPath="./RS_CutPng_Result/"#定义裁剪后的保存路径 forfileinfiles: a,b=os.path.splitext(file)#拆分影像图的文件名称 this_dir=os.path.join(dir+file) img=Image.open(this_dir)#按顺序打开某图片 width,hight=img.size w=480#宽度 h=360#高度 _id=1#裁剪结果保存文件名：0-N升序方式 y=0 while(y+h<=hight):#控制高度,图像多余固定尺寸总和部分不要了 x=0 while(x+w<=width):#控制宽度，图像多余固定尺寸总和部分不要了 new_img=img.crop((x,y,x+w,y+h)) new_img.save(ResultPath+a+"_"+str(_id)+b) _id+=1 x+=w y=y+h if__name__=='__main__': _path=r"./RS_TiffDir/"#遥感tiff影像所在路径 #裁剪影像图 Convert_To_Png_AndCut(_path)
将转换成png后的图加载到软件中（专业软件ENVI5.3）查看结果详细信息如下图所示，成功的转换成png格式了。
下面是加载裁剪后的影像图（Tiff格式的）
deftoCutPng(dir):函数效果图如下图所示。
以上这篇利用Python裁切tiff图像且读取tiff,shp文件的实例就是小编分享给大家的全部内容了，希望能给大家一个参考，也希望大家多多支持毛票票。
声明：本文内容来源于网络，版权归原作者所有，内容由互联网用户自发贡献自行上传，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任。如果您发现有涉嫌版权的内容，欢迎发送邮件至：czq8825#qq.com（发邮件时，请将#更换为@）进行举报，并提供相关证据，一经查实，本站将立刻删除涉嫌侵权内容。