python gdal Warp 矢量掩膜栅格

目录

矢量裁剪栅格代码

 cropToCutline 裁剪效果

 cropToCline 是否设置的优缺点

从矢量文件中选择部分要素，进行栅格裁剪生成与原始栅格大小相同的掩膜文件tif

其他参数对裁剪效果的影像

参考

---

矢量裁剪栅格代码

from osgeo import gdal,gdalconst
shppath = r'D:\Africa\Africa_city.shp'
tifpath = r'D:\regionImg\VNL_2012Africa.tif'
outtif1 = r'D:\Africa\Africa_FID0.tif'
cutlineWhere = 'FID = 2485'
ds = gdal.Warp(
outtif1,    #裁剪后图像保存的完整路径（包括文件名）
tifpath,    #待裁剪的影像完整路径（包括文件名）
format='GTiff', # 保存图像的格式
cutlineDSName=shppath, # 矢量文件的完整路径
cropToCutline=True, # 保证裁剪后影像大小跟矢量文件的图框大小一致（设置为False时，结果图像大小会跟待裁剪影像大小一样，则会出现大量的空值区域）
cutlineWhere=cutlineWhere #矢量文件筛选条件,
#dstNodata=0
)

 cropToCutline 裁剪效果

cropToCutline=True, 裁剪出的影像大小与矢量图像相近。  

边界像元若不完全在矢量图形内，则大概率被裁剪掉（这里还不确定边界像元什么条件下被裁剪掉，因为有的边界像元会被留下，有的则被裁剪掉。）。

关于掩膜像元（下图左，绿色部分），即裁剪后的栅格中不在矢量图形内的像元。若dstNodata不设置，则设为0，不作为Nodata。若dstNodata设置，则dstNodata设置的值作为Nodata的值，掩膜像元都设为nodata。dstNodata不设置，则裁剪后的栅格影像的nodata的设置取决于原始影像。此时若原始影像的nodata也没有设置，则掩膜像元值为0，否则掩膜像元为nodata，值为nodata对应的值。

cropToCutline=True (左：不设dstNodata，右：设dstNodata=-1)

dstNodata未设置 掩膜像元值为0

dstNodata=-1,掩膜像元为nodata

裁剪了两个tif

cropToCutline=False, 裁剪出的影像大小与原始栅格大小相同。不在矢量内的像元被掩膜。

cropToCutline=False (不设置dstNodata，绿色部分为掩膜像元)

裁剪后全局

缩放到矢量图像大小

FID=2473

 cropToCline 是否设置的优缺点

cropToCline=True	1）裁剪出的栅格较小，大小接近shp，裁剪速度较快，减小了后续加载、读、写、统计等的负担。 2）裁剪的栅格并不完全覆盖到矢量边界，与矢量边界之间有空隙。
cropToCline=False	1）裁剪出的栅格较大，与原始栅格相同。裁剪速度较慢。 2）裁剪的栅格也并不完全覆盖到矢量边界，有的图形裁剪会覆盖，有的没有覆盖（若上图右侧FID=2473这个图形裁剪出来的结果就没有覆盖。）

统计不同区域像素值，我会更倾向于用cropToCline=False。虽然它的计算量大，也不是所有矢量裁剪后都是全覆盖该矢量图形，但连个公共边矢量分别裁剪得到的栅格合并是无缝隙的，这样统计的结果可能更好些。而用cropToCline=True分别裁剪后进行统计，感觉统计结果会偏小一点。

从矢量文件中选择部分要素，进行栅格裁剪生成与原始栅格大小相同的掩膜文件tif

def getMaskTifByShp(shp_path,tifpath,outTifpath,sql):
    '''选择矢量文件中的部分要素，裁剪栅格，生成与输入栅格同等大小的mask.tif。
    生成结果中，像元值1为目标像元，像元值0为掩膜像元。'''
    #获取栅格信息
    inDs = gdal.Open(tifpath)
    rows = inDs.RasterYSize
    cols = inDs.RasterXSize
    geotrans = inDs.GetGeoTransform()
    proj = inDs.GetProjection()
    #创建内存栅格
    mem = gdal.GetDriverByName('MEM')
    mid_ds = mem.Create('', cols, rows, 1, gdal.GDT_Byte)
    mid_ds.GetRasterBand(1).WriteArray(np.ones((rows, cols), dtype=np.bool))
    mid_ds.SetGeoTransform(geotrans)
    mid_ds.SetProjection(proj)
    # #裁剪生成内存mask
    # mask_ds = gdal.Warp('', mid_ds, format='MEM', cutlineDSName=shp_path,cropToCutline=False,cutlineWhere=sql)
    # #输出
    # gtiff = gdal.GetDriverByName('GTiff')
    # result = gtiff.CreateCopy(outTifpath, mask_ds)
    # result.FlushCache()
    # del result,inDs,mid_ds,mask_ds
    #裁剪生成mask
    mask_ds = gdal.Warp(outTifpath, mid_ds, format='GTiff', cutlineDSName=shp_path,cropToCutline=False,cutlineWhere=sql)
    #输出
    del inDs,mid_ds,mask_ds


shp_path = r'D:\Africa_city.shp'
tifpath = r'D:\VNL_2012Africa.tif'
outTifpath = r'D:\2485.tif'
sql = 'FID = 2485'
getMaskTifByShp(shp_path,tifpath,outTifpath,sql)

 首先在内存中创建dataset，数据类型为byte。裁剪生成mask,根据format设置可以在内存中生成mask或者直接输出到硬盘。

输出结果：

白色为裁剪出的区域（白色DN=1，黑色为DN=0）

其他参数对裁剪效果的影响

1、固定cropToCline = True

1)指定裁剪后栅格的输出分辨率。

这里指定裁剪后的栅格分辨率与原始影像相同。在显示指定与不指定两种情况下，虽然输出影像的分辨率是相同的，但裁剪结果仍有差别（如下图）。

#指定xRes，yRes
mask_ds = gdal.Warp(outtifpath, tifpath, format='GTiff', cutlineDSName=shp_path, cropToCutline=True, cutlineWhere=sql,xRes=500,yRes=500)

两者叠加

未指定输出分辨率

指定输出分辨率

2）cutlineBlend  以像素为单位的剪切线混合距离。

mask_ds = gdal.Warp(outtifpath, tifpath, format='GTiff', cutlineDSName=shp_path, cropToCutline=True, cutlineWhere=sql,cutlineBlend=0.15,xRes=500,yRes=500)

cutlineBlend=None	cutlineBlend=0.1
cutlineBlend=0.5	cutlineBlend=0.2

 

参考

在内存中创建裁剪后的栅格参考自：Python3.GDAL从shp文件生成mask_碎积云-CSDN博客 

WarpOptions介绍：Python GDAL学习笔记(二)_江北20190411的博客-CSDN博客