【GDAL-Python】5-在Python中使用GDAL将栅格数据分割成相等的部分

1-介绍

1.1 主要内容

（1）在本教程中，使用GDAL将栅格数据拆分为几个相等的部分

（2）拆分手段： gdal.Warp或者 gdal.Translate

（3）关键流程：
1）找出较小瓦片的xmin、ymin、xmax和ymax坐标,（这里是以xmin等指的是地理坐标范围）
2） gdal.Warp或者 gdal.Translate进行裁剪

1.2 裁剪要点

（1）裁剪示意如下，主要找到每个小块的地理坐标范围xmin、ymin、xmax和ymax；
注意gt = dem.GetGeoTransform()中表示原始图像xmin = gt[0]，ymax = gt[3]，（ymin,xmin)是值的左上角坐标。

（2）gdal.Warp或者 gdal.Translate进行裁剪

• gdal.Warp(“dem”+str(i)+str(j)+“.tif”, dem, outputBounds = (xmin, ymin, xmax, ymax), dstNodata = -9999)
• gdal.Translate(“dem_translate”+str(i)+str(j)+“.tif”, dem, projWin = (xmin, ymax, xmax, ymin), xRes = res, yRes = -res)

2-代码实现

2.1 数据介绍

（1）原始数据信息：选用了1个波段的WordView3全色影像，影像分辨率为0.31米，影像元数据如下图：

（2）原始影像在QGIS显示效果：

2.2 代码实现

（1）计算裁剪分片
（2）进行裁剪

from osgeo import gdal

dem = gdal.Open(r"GDAL_testing_data\JAX_IMG1_PAN.TIF")
gt = dem.GetGeoTransform()

# ================================1.获取原始图像左上角的地理坐标===============================================
xmin = gt[0]+0*gt[1]+0*gt[2]#g[0-2]分别是x方向分辨率
ymin = gt[3]+0*gt[4]+0*gt[5]
res = gt[1]# 影像分辨率，注意x方向为正,y方向分辨率为负，但数值上x,y方向上一致


# ================================2.计算原始图像地理坐标范围长度和宽度========================================
xlen = res * dem.RasterXSize
ylen = res * dem.RasterYSize
print(f'xlen={xlen},ylen={ylen}')

# ================================3.设置x,y方向上的分块数=============================================
xdiv = 2
ydiv = 2


# ================================4.计算单块的尺寸====================================================
xsize = xlen/xdiv
ysize = ylen/ydiv
print(f'xsize={xsize},ysize={ysize}')

# ================================5.计算出分块的角点坐标===============================================
xsteps = [xmin + xsize * i for i in range(xdiv+1)]
print(f'xsteps={xsteps}')
ysteps = [ymin +ysize * i for i in range(ydiv+1)]
print(f'ysteps={ysteps}')

# ================================6.循环进行裁剪===============================================
for i in range(xdiv):
    for j in range(ydiv):
        xmin = xsteps[i]
        xmax = xsteps[i+1]
        ymin = ysteps[j]
        ymax = ysteps[j+1]
        
        print("xmin: "+str(xmin))
        print("xmax: "+str(xmax))
        print("ymin: "+str(ymin))
        print("ymax: "+str(ymax))
        print("\n")
        
        #使用 gdal warp裁剪
        # gdal.Warp("dem"+str(i)+str(j)+".tif", dem, 
        #           outputBounds = (xmin, ymin, xmax, ymax), dstNodata = -9999)
        #或者使用 gdal translate裁剪，但需要注意两个的传入裁剪范围顺序
        gdal.Translate("dem_translate"+str(i)+str(j)+".tif", dem, projWin = (ymin, xmin, ymax,xmax))
 
# 7.关闭
dem = None

2.3 均匀裁剪效果显示

改图亮度为人为显示，为了突出裁剪分块效果，左上和右下是经过2%线性拉伸显示，其他则是归一化后；
此外裁剪后会出现接边等问题，裁剪不充分等问题，有待解决（TODO)

3-参考资料

3.1 gdal.Translate两种裁剪方式

（1）一种是给定图像坐标系统下的坐标位置裁剪
（2）一种是给定地理坐标系统下的坐标位置裁剪

3.1.1 以给定图像坐标范围裁剪

图像坐标：x,y,w,h

执行代码参数：gdal.Translate(“dem_translate”+“.tif”, dem, srcWin = (y,x,h,w))，主要是srcWin表示传入图像坐标范围。
PS:注意其中srcWin以此传入y,x,h,w，一定不要传入x,y,w,h了，这里常有坑

3.1.2 以给定图像地理坐标范围裁剪

(1)图像转地理：为了验证，先把图像裁剪区域（x,y,w,h)转化为地理坐标，这里得理解影像得GeoTransform矩阵，转换函数如下

# 传入gdal.Open的gdal.Datast,以及图像坐标x,y
def GeoCoord(dataset, x, y):
    GT = dataset.GetGeoTransform()#(0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
    print(GT)
    if GT:
        return [GT[0] + x*GT[1] + y*GT[2], GT[3] + x*GT[4] + y*GT[5]]

（2）执行代码参数：gdal.Translate(“dem_translateProj”+“.tif”, dem, projWin = (ymin,xmin,ymax,xmax))，projWin表示传入地理坐标范围。
PS:这里的ymin,xmin,ymax,xmax分别用左下角坐标和右上角坐标计算，在这里表示如下

xmin, ymin=GeoCoord(dem,x,y)#计算图像左下角对应的地理坐标
print(f'xmin={xmin},ymin={ymin}')
xmax, ymax=GeoCoord(dem,x+w,y+h)# 计算图像右上角对应的地理坐标

3.1.3 代码实现

from osgeo import gdal
dem = gdal.Open(r"GDAL_testing_data\JAX_IMG1_PAN.TIF")
#=========================== 1.图像坐标裁剪区域=======================
x=1000
y=500
h=500
w=1000
# 按照图像坐标进行裁剪
gdal.Translate("dem_translateSrc"+".tif", dem, srcWin = (y,x,h,w))
#======================== 2.计算出图像坐标对应的地理坐标=====================
def GeoCoord(dataset, x, y):
    GT = dataset.GetGeoTransform()#(0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
    print(GT)
    if GT:
        return [GT[0] + x*GT[1] + y*GT[2], GT[3] + x*GT[4] + y*GT[5]]
# 2.1先计算对应的地理坐标
xmin, ymin=GeoCoord(dem,x,y)#xmin,ymin
print(f'xmin={xmin},ymin={ymin}')
xmax, ymax=GeoCoord(dem,x+w,y+h)
print(f'xmax={xmax},ymax={ymax}')
# 2.2 执行根据地理坐标范围进行裁剪
gdal.Translate("dem_translateProj"+".tif", dem, projWin = (ymin,xmin,ymax,xmax))#左下角和右上角

3.1.4 效果显示

影像裁剪范围为（x=1000，y=500，h=500，w=1000），从裁剪区域与整个原始图像的关系可以验证x＞y，w＞h。