GDAL：使用二进制栅格掩膜、栅格转矢量

使用二进制栅格掩膜

我们有代表数据掩码的二进制图像，我们可以将numpy.where（）函数与gdal_calc.py一起使用。假设我们有一个包含像素值0和1的binary.tif。我们现在想将binary.tif=0时，对应位置的merged.tif中的所有像素设置为无数据，我们想使用-32768作为NoData值，可以使用如下命令：

gdal_calc.py -A merged.tif -B binary.tif --calc="numpy.where(B==1, A, -32768)"  --NoDataValue -32768 --outfile masked.tif 

该命令指出，只要B栅格的值为1，输出就应该是来自A的像素值，否则应该是-32768。因为我们还将-NoDataValue设置为-32768，所以该命令有效地设置了条件不匹配的所有像素(B！= 1)到NoData。这里对应的就是ArcGIS栅格计算器的SetNull函数。

栅格转矢量

gdal_polygonize.py支持栅格转矢量，假设我们从之前案例中，合并切片部分的合并栅格中提取最高高程的坐标。使用gdalinfo -stats查询栅格显示最高像素值为8748，我们可以使用gdal_calc.py来创建一个栅格，使用的条件是只匹配具有该值的像素。

cd D:\GDALLearn\gdal-tools\srtm
python C:\Users\YIN\miniconda\Scripts\gdal_calc.py --calc 'A==8748' -A merged.vrt --outfile everest.tif --NoDataValue=0

上面这样的布尔表达式的结果将是一个像素值为1和0的光栅。当我们将 NoData 设置为0时，条件匹配的值为1的像素只有1个，再可以使用 gdal _ polygonize.py 将其转换为矢量shp文件。

python C:\Users\YIN\miniconda\Scripts\gdal_polygonize.py everest.tif everest.shp

如果要提取多边形的质心并打印坐标，可以使用ogrinfo命令。

ogrinfo everest.shp -sql 'SELECT AsText(ST_Centroid(geometry)) from everest' -dialect SQLite