| python gdal Warp 矢量掩膜栅格 | 您所在的位置：网站首页 › gdal生成tiff图例 › python gdal Warp 矢量掩膜栅格 |
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=False1）裁剪出的栅格较大，与原始栅格相同。裁剪速度较慢。 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博客 |
【本文地址】