前几天在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写代码跑的话，效率就低很多，可能还有改进的地方吧。