前幾天在DMSP數據預處理中介紹了利用gdal.warp()來實現影像裁剪，今天介紹如何通過代碼一步步實現影像裁剪功能。影像裁剪的過程其實就是讀取影像數據，獲取影像數據的仿射、投影信息，將地理坐標轉換成像素坐標，讀取用于裁剪的矢量數據，獲取相應的投影信息，將矢量數據范圍坐標轉換成像素坐標，生成用于裁剪的掩膜文件(0，1)，根據掩膜文件，生成掩膜范圍內實際裁剪array(可以理解為矩陣點乘，掩膜文件為1的保留有影像數據值，為0的則影像數據值變為0)，最后對數據進行寫入，生成裁剪后的影像文件。

結合實際操作，逐步說明吧。

首先第一步影像數據、柵格數據的讀入，及相關信息的生成(主要就是放射信息，地理坐標轉換為像素坐標)。先是影像數據讀入。

#定義基礎數據打開及基礎信息獲取函數

def Get_Imgdata(filepath):

ds = gdal.Open(filepath)

#img_band = ds.GetRasterBand(1)

#band_arr = img_band.ReadAsArray()

#bandnum = img_ds.RasterCounter

#img_width =img_band.XSize

#img_height = img_band.YSize

#img_bandnum = img_ds.RasterCount

geo = ds.GetGeoTransform()

prj = ds.GetProjection()

return ds,geo,prj

影像數據讀入及相關基礎信息獲取還是比較容易的，接下來是矢量數據讀入和坐標轉換。

def Get_shpdata(filepath,img_geo):

ds = ogr.Open(filepath,0)

layer = ds.GetLayer()

shp_extent = layer.GetExtent()

#luxy = (shp_extent[0],shp_extent[3])

ulx,uly = shp_extent[0],shp_extent[3]

#rlxy = (shp_extent[1],shp_extent[2])

lrx,lry = shp_extent[1],shp_extent[2]

#實際坐標轉換成像素坐標，獲取像素裁剪范圍

inv_img_geo = gdal.InvGeoTransform(img_geo)

ulxy_pixel = World2Pixel(inv_img_geo,ulx,uly)

lrxy_pixel = World2Pixel(inv_img_geo,lrx,lry)

#pxwidth = lrxy_pixel[0]-ulxy_pixel[0]

#pxheight = lrxy_pixel[1]-ulxy_pixel[1]

print(ulxy_pixel,lrxy_pixel)

feat = layer.GetFeature(0)

feat_geom = feat.GetGeometryRef()

points = feat_geom.GetGeometryRef(0)

points_list = []

points_pixel =[]

for p in range(points.GetPointCount()):

points_list.append((points.GetX(p),points.GetY(p)))

#轉換成像素坐標

#在轉之前需要修改仿射信息，不修改的話后面的mask文件會出錯

img_geo0 = list(img_geo)

img_geo0[0] = ulx

img_geo0[3] = uly

inv_img_geo0 = gdal.InvGeoTransform(img_geo0)

for p in points_list:

points_pixel.append(World2Pixel(inv_img_geo0,p[0],p[1]))

return ulxy_pixel,lrxy_pixel,img_geo0,points_pixel

def World2Pixel(geotrans,x,y):

#inv_geotrans = gdal.InvGeoTransform(geotrans)

xy_pixel = gdal.ApplyGeoTransform(geotrans,x,y)

# 將轉換后的結果進行取整，必須要做的

int_x_pixel,int_y_pixel = map(int,xy_pixel)

return (int_x_pixel,int_y_pixel)

這一部分內容就比較復雜一些，里面有一些小的技巧和一些容易出錯的地方。第一個要說明的是我定義了一個地理坐標轉換像素坐標的函數，在這里很多其他類似的操作，都是利用仿射信息進行計算(仿射信息內提供了數據的起始點，也就是左上角的坐標，根據pixelwidth，pixelheight進行計算)，而我給大家提供更為方便的兩個函數，gdal.InvGeoTransform和ApplyGeoTransform，可以輕松的實現地理坐標與像素坐標之間的轉換。第二個要說明的是矢量拐點坐標的獲取，對于規則形狀的矢量范圍來說，很容易，但是實際情況是，裁剪矢量范圍一般是不規則的，比如一個地方的行政區劃，這樣的矢量數據通常含有眾多拐點，因此獲取每一個拐點坐標，并進行坐標轉換，才能實現后面的裁剪工作。 feat = GetGeometryRef() 和points = feat_geom.GetGeometryRef(0)結合循環，就可以逐個提取矢量要素的拐點坐標了。第三個要說的是仿射信息的更改。因為是需要用矢量數據進行裁剪，后續的掩膜文件生成和文件的寫入，的仿射信息的起始點坐標應該是矢量文件的左上角坐標，這是需要修改的，也是容易出錯的地方。

再接下來就是掩膜文件的生成。生成一個范圍與矢量數據一致的淹沒文件，掩膜文件范圍是矢量數據的四至范圍，矢量數據范圍內的掩膜數值為1，之外的為0，大家可以去自己生成試一下。

def imageToArray(i):

a=gdalnumeric.fromstring(i.tobytes(),'b')

a.shape=i.im.size[1], i.im.size[0]

return a

def Create_Mask(x,y,p):

clippoly = Image.new("L", (x,y), 1)

draw_ploy = ImageDraw.Draw(clippoly)

draw_ploy.polygon(p, 0)

mask = imageToArray(clippoly)

print(mask)

return mask

然后就是實現數據的裁剪。

def Clip(choose,ds,mask,new_geo,columns,rows,offset):

band = ds.GetRasterBand(1)

driver = gdal.GetDriverByName('GTiff')

out_ds = driver.Create(r'D:\gdal_data\CLIP\B4_Clip.tif',columns,rows,1,band.DataType)

out_ds.SetProjection(img_ds.GetProjection())

out_ds.SetGeoTransform(new_geo)

out_band = out_ds.GetRasterBand(1)

band_arr = band.ReadAsArray(offset[0],offset[1],columns,rows)

if choose == 'v':

clip_data = gdalnumeric.choose(mask, (band_arr,0))

if choose == 'c':

clip_data = band_arr

out_band.WriteArray(clip_data)

out_ds.FlushCache()

del out_ds

我這里寫著玩的，'V'代表不規則矢量數據，'C'代表規則矢量數據。體驗一下。最后就是調用前面所有操作的主函數了。

if __name__ == '__main__':

img_path = r'D:\gdal_data\CLIP\LC81270372014355LGN00_B4.tif'

shp_path = r'D:\gdal_data\CLIP\范圍_proj.shp'

img_ds,img_geo,img_prj = Get_Imgdata(img_path)

UL_pixel,LR_pixel,img_geonew,pts_pixel = Get_shpdata(shp_path,img_geo)

cols = LR_pixel[0] - UL_pixel[0]

rows = LR_pixel[1] - UL_pixel[1]

clip_mask = Create_Mask(cols,rows,pts_pixel)

clip_choose = input()

if clip_choose == 'v':

print('Execute Vector Clip!!')

Clip(clip_choose,img_ds,clip_mask,img_geonew,cols,rows,UL_pixel)

if clip_choose == 'c':

print('Execute Regular_shape Clip!!')

Clip(clip_choose,img_ds,clip_mask,img_geonew,cols,rows,UL_pixel)

else:

print("please enter 'v' or 'c',otherwise there is no way to clip!" )

#print()

#Clip()

至此，就實現了數據的裁剪。照例貼上裁剪的案例吧。這里我只裁剪了單一波段，起始多波段裁剪過程一樣，就是最后寫入的時候變換一下數據寫入的方式，前面介紹過了，就不再重復說了。

OK，今天就分享這些。之前說過要簡單介紹一下GEE的python的一些基礎，例如ee.Geometry，ee.Feature, ee.Image等，我想在介紹這個內容之前，還想簡單介紹一下gdal矢量數據創建，就比如如何創建點、線、面等，之后再介紹GEE，這樣能夠有個比較好的對比。其實GEE的話，還是推薦大家用google自帶的，本地python化的話，有個很大的問題就是效率的問題。我用GEE平臺做一個省的土地利用分類工作，在基礎數據準備好的基礎上，用GEE自帶的分類模型跑的話，非常的快，而我用python寫代碼跑的話，效率就低很多，可能還有改進的地方吧。

總結

以上是生活随笔為你收集整理的pythonmapdel_地质男转行学遥感Python——遥感数据裁剪的具体实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。